Housekeeping

Preamble

By Katherine Duchesneau

We are working with the roots Lesion and Mycorrhiza dataset for now.

R details

This section contains all the necessary details about the R version and package version used to run this Notebook.

  • File creation date: 2018-02-12
  • R version 3.4.1 (2017-06-30)
  • tidyverse package version: 1.2.1
  • lme4 package version 1.1.13
  • MASS package version 7.3.47
  • MuMIn package version 1.15.6
  • emmeans package version 0.9.1
  • boot package version 1.3.20
  • brms package version 2.0.1
  • loo package version 1.1.0
  • fitdistrplus package version 1.0.9
  • mefa package version 3.2.7
  • library(devtools) and install_github("vqv/ggbiplot") package version 0.55
  • vegan package version 2.4.3
  • glmmADMB package version 0.8.3.3
  • car package version 2.1.6
  • factoextra package version 1.0.5.999
  • pvclust package version 2.0.0

Load Packages and Functions

library(tidyverse)
library(lme4)
library(MASS)
library(MuMIn)
library(emmeans)
library(boot)
library(brms)
library(loo)
library(fitdistrplus)
library(mefa)
library(ggbiplot)
library(pvclust)
library(vegan)
library(ggbiplot)
library(factoextra)
library(glmmADMB)
library(car)

Load a custome theme for clean, readable graphs:

Overdispersion function, tests for overdispersion in GLMM:

Introducing the Data

Load/ Modify data: Co-occurence

## [1] 99.55752
## [1] 100
##      indiv2       None.Path           None.Myc           Herbivory      
##  10IAH1 :  100   Length:22500       Length:22500       Min.   :0.00000  
##  10IAH2 :  100   Class :character   Class :character   1st Qu.:0.00000  
##  10ISC1 :  100   Mode  :character   Mode  :character   Median :0.00000  
##  10ISC2 :  100                                         Mean   :0.03658  
##  10ISC3 :  100                                         3rd Qu.:0.00000  
##  10OAH1 :  100                                         Max.   :1.00000  
##  (Other):21900                                                          
##    None.Herb          Lesion        Mycorrhiza    
##  Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:1.0000   1st Qu.:0.000   1st Qu.:0.0000  
##  Median :1.0000   Median :1.000   Median :1.0000  
##  Mean   :0.9631   Mean   :0.622   Mean   :0.5732  
##  3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.000   Max.   :1.0000  
## 
##     indiv2            None.Path        None.Myc        Herbivory      
##  Length:22500       Min.   :0.000   Min.   :0.0000   Min.   :0.00000  
##  Class :character   1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Mode  :character   Median :0.000   Median :0.0000   Median :0.00000  
##                     Mean   :0.378   Mean   :0.4268   Mean   :0.03658  
##                     3rd Qu.:1.000   3rd Qu.:1.0000   3rd Qu.:0.00000  
##                     Max.   :1.000   Max.   :1.0000   Max.   :1.00000  
##                                                                       
##    None.Herb          Lesion        Mycorrhiza          pop      
##  Min.   :0.0000   Min.   :0.000   Min.   :0.0000   5      :2900  
##  1st Qu.:1.0000   1st Qu.:0.000   1st Qu.:0.0000   7      :2600  
##  Median :1.0000   Median :1.000   Median :1.0000   8      :2600  
##  Mean   :0.9634   Mean   :0.622   Mean   :0.5732   3      :2200  
##  3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:1.0000   4      :2200  
##  Max.   :1.0000   Max.   :1.000   Max.   :1.0000   1      :2000  
##                                                    (Other):8000  
##  location  species  
##  I:10600   AH:3100  
##  O:11900   CL:3900  
##            GA:4900  
##            GR:1300  
##            MR:3100  
##            SC:6200  
## 
##     indiv2          None.Path None.Myc  Herbivory None.Herb Lesion   
##  Length:22500       0:13995   0:12897   0:21677   0:  823   0: 8505  
##  Class :character   1: 8505   1: 9603   1:  823   1:21677   1:13995  
##  Mode  :character                                                    
##                                                                      
##                                                                      
##                                                                      
##                                                                      
##  Mycorrhiza      pop       location  species  
##  0: 9603    5      :2900   I:10600   AH:3100  
##  1:12897    7      :2600   O:11900   CL:3900  
##             8      :2600             GA:4900  
##             3      :2200             GR:1300  
##             4      :2200             MR:3100  
##             1      :2000             SC:6200  
##             (Other):8000
## [1] 225

Load/ Modify data: Totals

## [1] 226
## [1] 226
##      indiv         Decay           Pathogen          Hyphae     
##  10IAH1 :  1   Min.   : 0.000   Min.   : 0.000   Min.   : 2.00  
##  10IAH2 :  1   1st Qu.: 1.000   1st Qu.: 2.000   1st Qu.:34.00  
##  10ISC1 :  1   Median : 4.000   Median : 4.000   Median :48.00  
##  10ISC2 :  1   Mean   : 7.894   Mean   : 7.168   Mean   :47.23  
##  10ISC3 :  1   3rd Qu.: 9.000   3rd Qu.:10.000   3rd Qu.:63.00  
##  10OAH1 :  1   Max.   :60.000   Max.   :86.000   Max.   :96.00  
##  (Other):220                                                    
##    None.Path       Arbuscules        Vesicles         M_Hyphae    
##  Min.   : 0.00   Min.   : 0.000   Min.   : 0.000   Min.   : 0.00  
##  1st Qu.:20.00   1st Qu.: 0.000   1st Qu.: 1.000   1st Qu.:30.25  
##  Median :33.00   Median : 0.000   Median : 6.000   Median :46.00  
##  Mean   :37.71   Mean   : 2.363   Mean   : 9.956   Mean   :45.16  
##  3rd Qu.:53.75   3rd Qu.: 3.000   3rd Qu.:14.000   3rd Qu.:61.00  
##  Max.   :97.00   Max.   :21.000   Max.   :73.000   Max.   :94.00  
##                                                                   
##     None.Myc        Herbivory        None.Herb           pop     location
##  Min.   :  0.00   Min.   : 0.000   Min.   : 62.00   5      :29   I:106   
##  1st Qu.: 18.00   1st Qu.: 0.000   1st Qu.: 96.00   8      :27   O:120   
##  Median : 38.50   Median : 1.000   Median : 99.00   7      :26           
##  Mean   : 42.52   Mean   : 3.624   Mean   : 96.38   3      :22           
##  3rd Qu.: 66.00   3rd Qu.: 4.000   3rd Qu.:100.00   4      :22           
##  Max.   :100.00   Max.   :38.000   Max.   :100.00   1      :20           
##                                                     (Other):80           
##  species
##  AH:31  
##  CL:39  
##  GA:50  
##  GR:13  
##  MR:31  
##  SC:62  
## 

Field survey data

FloristicSurvey <-read.csv("CSV/May_2017_FloristicSurvey_summer2016.csv", sep="\t")
PlantLength <- read.csv("CSV/May_2017_PlantLength_summer2016.csv", sep="\t")
Soil_characteristics <- read.csv("CSV/Soil_characteristics.txt", sep="\t")

#Restructuring
head(PlantLength)
str(PlantLength)
## 'data.frame':    231 obs. of  4 variables:
##  $ Pop         : int  14 14 14 14 14 14 14 14 14 14 ...
##  $ Location    : Factor w/ 2 levels "in","out": 2 2 2 2 2 1 1 1 1 1 ...
##  $ Species     : Factor w/ 7 levels "Acer_saccharum",..: 1 4 3 5 7 1 4 3 5 7 ...
##  $ Length_plant: num  14 87 20 12 67 11 65 51 20 100 ...
colnames(PlantLength)[1] <- "Population"
PlantLength$Population<-as.factor(PlantLength$Population)
PlantLength$Location<-as.factor(PlantLength$Location)

FloristicSurvey$Population<-as.factor(FloristicSurvey$Population) # Make the population a factor and not an interger
FloristicSurvey[is.na(FloristicSurvey)] <- 0
sum(row.has.na <- apply(FloristicSurvey, 1, function(x){any(is.na(x))}))
## [1] 0

Root size

Root_size <-read.csv("CSV/Root_size.csv")
Root_size$pop<-gsub("^([0-9]+)[I,O][A-Z]+.*","\\1",Root_size$indiv)
Root_size$location<-gsub("^[0-9]+([I,O])[A-Z]+.*","\\1",Root_size$indiv)
Root_size$species<-gsub("^[0-9]+[I,O]([A-Z]+).*","\\1",Root_size$indiv)
Root_size$location<-as.factor(Root_size$location)
Root_size$species<-as.factor(Root_size$species)
Root_size$pop<-as.factor(Root_size$pop)

Variable descriptions

The variables are:

Co-occurence:

  • indiv2: The unique code given to a sample. In this dataset the unique sample code repeats 100 times for each cross where an observation was taken on the root sample.

  • None.Path: A binary representation of whether a sign of pathogen activity was recorded (0) or not (1).

  • Lesion: The counterpart to the previous variable (None.Path). A binary representation of whether a sign of pathogen activity was recorded (1) or not (0).

  • None.Myc: A binary representation of whether a sign of myccorhizal activity was recorded (0) or not (1).

  • Mycorrhiza: The counterpart to the previous variable (None.Myc). A binary representation of whether a sign of mycorrhizal activity was recorded (1) or not (0).

  • None.Herb: A binary representation of whether herbivory was recorded (0) or not (1).

  • Herbivory: The counterpart to the previous variable (None.Herb). A binary representation of whether Herbivoryn was recorded (1) or not (0).

  • Population (pop): The coding number representing the population at which the sample was collected.

  • location: A code representing the whether the sample was collected inside a Alliaria petiolata population (I) or whther it was collected at least 7 m outside of he furthest individual in the A. petiolata population (O).

  • species: The particular species to which the sample belongs.

  • Cross: The particular cross number where the data was recorded on the individual sample.

Total scores:

  • indiv: The unique code given to a sample.

  • Decay: The total number of crosses where decay was recorded an individual sample when doing a total of 100 crosses/ sample.

  • Pathogen: The total number of crosses where pathogen was recorded on an individual sample when doing a total of 100 crosses/ sample.

  • Hyphae: The total number of crosses where non-myccorhizal hyphae was recorded on an individual sample when doing a total of 100 crosses/ sample.

  • None.Path: The total number of crosses where no signs of pathogen activities were recorded on an individual sample when doing a total of 100 crosses/ sample. Note that the total of Decay, Pathogen, Hyphae, and None.Path must come to 100 per individual to account for all 100 crosses.

  • Arbuscule: The total number of crosses where an arbuscule was recorded on an individual sample when doing a total of 100 crosses/ sample.

  • Vesicules: The total number of crosses where a vesicule was recorded on an individual sample when doing a total of 100 crosses/ sample.

  • M_Hyphae: The total number of crosses where myccorhizal hyphae was recorded on an individual sample when doing a total of 100 crosses/ sample.

  • None.Myc: The total number of crosses where no signs of mycorrhizal activities were recorded on an individual sample when doing a total of 100 crosses/ sample. Note that the total of Arbuscules, Vesicules, M_Hyphae, and None.Myc must come to 100 per individual to account for all 100 crosses.

  • Herbivory: The total number of crosses where herbivory was recorded on an individual sample when doing a total of 100 crosses/ sample.

  • None.Herb: The total number of crosses where no signs of herbivory was recorded on an individual sample when doing a total of 100 crosses/ sample. Note that the total of Herbivory, and None.Herb must come to 100 per individual to account for all 100 crosses.

  • Population (pop): The coding number representing the population at which the sample was collected.

  • location: A code representing the whether the sample was collected inside a Alliaria petiolata population (I) or whther it was collected at least 7 m outside of he furthest individual in the A. petiolata population (O).

  • species: The particular species to which the sample belongs.

Reviewing of hypotheses:

Question 1: Is there evidence that Mycorrhizal colonization provides protection against root pathogen damage?

Question 2: Do these effects vary among species?

Hypothesis 1

Mycorrhizal colonization protects roots from lesions and herbivory.

Hypothesis 2

A. petiolata invasion reduces mycorrhizal colonization and increases pathogen colonization and herbivory

Hypothesis 3

Cycorrhizal colonization is reduced in the presence of A. petiolata

Hypothesis 4

Mycorrhizal protection is completely lost in the presence of A. petiolata

Molecular data Hypothesis

Hypothesis 1

A. petiolata invasion changes the diversity of mycorrhizal and other fungal species present in the soil.

I will find that the molecular data from the soil extraction supports this claim by having a reduced biodiversity and abundance of mycorrhiza in samples taken inside A. petiolata populations.

Hypothesis 2

There will be a negative correlation between mycorrhizal and fungal pathogen diversity.

Field Data Analysis

Plant length and GM presence

row.has.na <- apply(PlantLength, 1, function(x){any(is.na(x))})
sum(row.has.na) 
## [1] 2
PlantLength.filtered <- PlantLength[!row.has.na,]# Removed all the rows with NAs
PlantLength.filtered$Pop <- as.factor(PlantLength.filtered$Pop)

car::qqp(PlantLength$Length_plant, "norm")

shapiro.test(PlantLength$Length_plant)
## 
##  Shapiro-Wilk normality test
## 
## data:  PlantLength$Length_plant
## W = 0.89321, p-value = 1.14e-11
mod1<-lm(Length_plant~Location*Species, data=PlantLength.filtered)

bc<-boxcox(mod1, lambda = seq(-2, 2, 0.1))

lambda <- bc$x[which.max(bc$y)]
PlantLength.filtered <- cbind(PlantLength.filtered, ((PlantLength.filtered$Length_plant^lambda)-1)/lambda)
names(PlantLength.filtered)[length(PlantLength.filtered)] <- "Yprime"
mod1BC <- lm(Yprime~Location*Species, data=PlantLength.filtered) 

par(mfrow = c(2, 2))
plot(mod1BC)

# This is basically a log transformation. 

# Multicolinearity?
car::vif(mod1BC)
##                       GVIF Df GVIF^(1/(2*Df))
## Location           5.43980  1        2.332338
## Species           94.60756  6        1.461035
## Location:Species 306.05628  6        1.611202
sqrt(car::vif(mod1BC)) > 2
##                  GVIF    Df GVIF^(1/(2*Df))
## Location         TRUE FALSE           FALSE
## Species          TRUE  TRUE           FALSE
## Location:Species TRUE  TRUE           FALSE
mod1<-lmer(log(Length_plant)~Location*Species+(1|Pop), data=PlantLength.filtered)

mod2<-lmer(log(Length_plant)~Location+(1|Pop), data=PlantLength.filtered)

mod3<-lmer(log(Length_plant)~1+(1|Pop), data=PlantLength.filtered)

AICres <- model.sel(mod1,mod2,mod3)
r.squaredLR(mod1)
## [1] 0.7640355
## attr(,"adj.r.squared")
## [1] 0.8900227
summary(mod1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(Length_plant) ~ Location * Species + (1 | Pop)
##    Data: PlantLength.filtered
## 
## REML criterion at convergence: 161.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.94568 -0.62634 -0.05365  0.56828  2.80893 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Pop      (Intercept) 0.006596 0.08122 
##  Residual             0.100384 0.31683 
## Number of obs: 229, groups:  Pop, 7
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                               2.34167    0.07565  30.955
## Locationout                              -0.03789    0.09778  -0.388
## SpeciesAnemone_hepatica                   0.23682    0.12794   1.851
## SpeciesCircaea_lutetiana                  1.08300    0.10771  10.055
## SpeciesGallium_aparine                    1.46347    0.10204  14.342
## SpeciesGeranium_robertinum                0.80086    0.10771   7.435
## SpeciesMaianthemum_racemosum              0.55442    0.12295   4.509
## SpeciesSolidago_canadensis                1.63655    0.09778  16.738
## Locationout:SpeciesAnemone_hepatica       0.21256    0.17101   1.243
## Locationout:SpeciesCircaea_lutetiana     -0.24113    0.14797  -1.630
## Locationout:SpeciesGallium_aparine       -0.20452    0.14392  -1.421
## Locationout:SpeciesGeranium_robertinum    0.07377    0.14799   0.498
## Locationout:SpeciesMaianthemum_racemosum  0.12277    0.16743   0.733
## Locationout:SpeciesSolidago_canadensis   -0.11605    0.13828  -0.839
## 
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
confint(mod1)
## Computing profile confidence intervals ...
##                                                2.5 %     97.5 %
## .sig01                                    0.01963921 0.16258131
## .sigma                                    0.28096327 0.33850558
## (Intercept)                               2.19817644 2.48515490
## Locationout                              -0.22469501 0.14891543
## SpeciesAnemone_hepatica                  -0.01319870 0.47949790
## SpeciesCircaea_lutetiana                  0.87828687 1.29033635
## SpeciesGallium_aparine                    1.26920608 1.65925061
## SpeciesGeranium_robertinum                0.59591108 1.00763645
## SpeciesMaianthemum_racemosum              0.31810708 0.78817334
## SpeciesSolidago_canadensis                1.44974225 1.82335269
## Locationout:SpeciesAnemone_hepatica      -0.11265649 0.54145529
## Locationout:SpeciesCircaea_lutetiana     -0.52495034 0.04066983
## Locationout:SpeciesGallium_aparine       -0.47948584 0.07045383
## Locationout:SpeciesGeranium_robertinum   -0.20959317 0.35592002
## Locationout:SpeciesMaianthemum_racemosum -0.19718336 0.44252670
## Locationout:SpeciesSolidago_canadensis   -0.38023461 0.14813034

Floristic composition and GM presence

Linear discriminate analysis

# Discriminant function analysis is used (general LDA format) because i have a cat. dependent and a binary independent. (it needs binary or continuous!)
# you might have to bin covers because it only works with categorical, so 0-5% none, 6-40 low etc.
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,25,50,75,100), labels = c("None", "Med_Low", "Med_High", "Overtaken"))
table(FloristicSurvey$GM_Coverage_category)
## 
##      None   Med_Low  Med_High Overtaken 
##        82        24        58        34
# [0-20] None
# ]20-50] Medium low
# ]50-80] Medium high
# ]80-100] Overtaken
#LDA
FloristicSurveysmall<-data.frame(FloristicSurvey[,4:56])
# Linear discriminant function analysis

FloristicSurveysmall <- FloristicSurveysmall[,-c(42,45)] 
FloristicSurveyLDA <- lda(GM_Coverage_category ~., data=FloristicSurveysmall)
# Extract scaling vectors
scalvec<-data.frame(FloristicSurveyLDA$scaling)
# Extract predictions
FloristicSurveyLDAval <- data.frame(predict(FloristicSurveyLDA)$x)

# Plot results
ggplot(data=FloristicSurveyLDAval,aes(x=LD1, y=LD2, group=FloristicSurveysmall$GM_Coverage_category))+
  stat_ellipse(geom="polygon",aes(colour=FloristicSurveysmall$GM_Coverage_category),fill=NA,size=1.2,alpha=0.3)+
  stat_ellipse(geom="polygon",aes(fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=1.2,alpha=0.3)+
  geom_point(aes(shape=FloristicSurveysmall$GM_Coverage_category,fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=I(4),alpha=I(0.8))+
  xlab("LD Axis 1")+ylab("LD Axis 2")+theme_simple() +
  theme(legend.title=element_blank())

#test significance of axis
anova(lm(FloristicSurveyLDAval$LD1~FloristicSurveysmall$GM_Coverage_category))
anova(lm(FloristicSurveyLDAval$LD2~FloristicSurveysmall$GM_Coverage_category))
#LDA WITHOUT A. PETIOLATA
FloristicSurveysmall<-data.frame(FloristicSurvey[,4:56])
# Linear discriminant function analysis
FloristicSurveysmall[10]<-NULL
FloristicSurveysmall[c(41,44)] <- NULL 
FloristicSurveyLDA<-lda(GM_Coverage_category ~ ., data=FloristicSurveysmall)
# Extract scaling vectors
scalvec<-data.frame(FloristicSurveyLDA$scaling)
# Extract predictions
FloristicSurveyLDAval <- data.frame(predict(FloristicSurveyLDA)$x)

# Plot results
ggplot(data=FloristicSurveyLDAval,aes(x=LD1,y=LD2,group=FloristicSurveysmall$GM_Coverage_category))+
  stat_ellipse(geom="polygon",aes(colour=FloristicSurveysmall$GM_Coverage_category),fill=NA,size=1.2,alpha=0.3)+
  stat_ellipse(geom="polygon",aes(fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=1.2,alpha=0.3)+
  geom_point(aes(shape=FloristicSurveysmall$GM_Coverage_category,fill=FloristicSurveysmall$GM_Coverage_category,colour=FloristicSurveysmall$GM_Coverage_category),size=I(4),alpha=I(0.8))+
  xlab("LD Axis 1")+ylab("LD Axis 2")+theme_simple() +
  theme(legend.title=element_blank())

#test significance of axis
anova(lm(FloristicSurveyLDAval$LD1~FloristicSurveysmall$GM_Coverage_category))
anova(lm(FloristicSurveyLDAval$LD2~FloristicSurveysmall$GM_Coverage_category))

Cluster analysis

# Ward Hierarchical Clustering with Bootstrapped p values
#do they cluster by population?
FloristicSurvey[is.na(FloristicSurvey)] <- 0
mydata<-t(data.frame(FloristicSurvey[,1],FloristicSurvey[,4:55]))
my.names <- mydata[1,]
colnames(mydata) <- my.names
mydata <- mydata[-1,]
d <- dist(t(mydata), method = "euclidean") # distance matrix
fit <- hclust(d, method = "ward.D")
fitd <- as.dendrogram(fit)

par(mfrow=c(3,2))
plot(fit,cex=0.4,main="Clustering with populations")
plot(cut(fitd, h=10)$upper, main="Upper tree of cut at h=10")
plot(cut(fitd, h=10)$lower[[1]], main="First branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[2]], main="Second branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[3]], main="Third branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[4]], main="Fourth branch of lower tree with cut at h=10")

#do they cluster by GM coverage?
FloristicSurvey[is.na(FloristicSurvey)] <- 0
#this will just make the dendogram easier to see
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,25,50,75,100), labels = c("None", "ML", "MH", "High"))
mydata<-t(data.frame(FloristicSurvey[,56],FloristicSurvey[,4:55]))
my.names <- mydata[1,]
colnames(mydata) <- my.names
mydata <- mydata[-1,]
d <- dist(t(mydata), method = "euclidean") # distance matrix
fit <- hclust(d, method = "ward.D")
fitd <- as.dendrogram(fit)
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,5,30,50,70,95,100), labels = c("None", "Low", "Med_Low", "Med_High", "High", "Overtaken"))

par(mfrow=c(3,2))
plot(fit,cex=0.4,main="Clustering with populations")
plot(cut(fitd, h=10)$upper, main="Upper tree of cut at h=10")
plot(cut(fitd, h=10)$lower[[1]], main="First branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[2]], main="Second branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[3]], main="Third branch of lower tree with cut at h=10")
plot(cut(fitd, h=10)$lower[[4]], main="Fourth branch of lower tree with cut at h=10")

NMDS

# https://jonlefcheck.net/2012/10/24/nmds-tutorial-in-r/
# https://oliviarata.wordpress.com/2014/04/17/ordinations-in-ggplot2/
#Make a matrix with no row or column equal to 0 (do not enclude the env variable (GM COVERAGE))
FloristicSurvey$GM_Coverage_category<-cut(FloristicSurvey$GM_Coverage, c(-Inf,25,50,75,100), labels = c("Low", "Med_Low", "Med_High", "High"))
table(FloristicSurvey$GM_Coverage_category)
## 
##      Low  Med_Low Med_High     High 
##       82       24       58       34
M <- as.matrix(FloristicSurvey[1:198,4:56])
M[is.na(M)] <- 0
rownames(M) <- FloristicSurvey$Population
M<-M[-c(130,133),]
M<-M[,-45]
GM_coverage_df <- data.frame(M[,52])
M<-M[,-52]
class(M) <- "numeric"


# with vegdist from Bray to seroeson: add binary = T 
dist_FloristicSurvey <- vegdist(M, method = "bray", binary = T)
#The metaMDS analysis could have done the distance matrix internally but i would rather control it since i have presence/abscence
meta.nmds.FloristicSurvey2D <- metaMDS(dist_FloristicSurvey, k=2, trymax = 1000)
## Run 0 stress 0.2097898 
## Run 1 stress 0.2103351 
## Run 2 stress 0.2099173 
## ... Procrustes: rmse 0.01192537  max resid 0.1277007 
## Run 3 stress 0.2126692 
## Run 4 stress 0.2131916 
## Run 5 stress 0.2130193 
## Run 6 stress 0.2186913 
## Run 7 stress 0.2098392 
## ... Procrustes: rmse 0.009869736  max resid 0.1035989 
## Run 8 stress 0.2352097 
## Run 9 stress 0.211495 
## Run 10 stress 0.2097224 
## ... New best solution
## ... Procrustes: rmse 0.006331696  max resid 0.0529704 
## Run 11 stress 0.2103191 
## Run 12 stress 0.2122882 
## Run 13 stress 0.2101326 
## ... Procrustes: rmse 0.009866801  max resid 0.09557844 
## Run 14 stress 0.2099405 
## ... Procrustes: rmse 0.01634622  max resid 0.1525454 
## Run 15 stress 0.2129822 
## Run 16 stress 0.2102758 
## Run 17 stress 0.2133365 
## Run 18 stress 0.2102033 
## ... Procrustes: rmse 0.006775596  max resid 0.07193472 
## Run 19 stress 0.2097459 
## ... Procrustes: rmse 0.008284148  max resid 0.09607869 
## Run 20 stress 0.2117128 
## Run 21 stress 0.2131278 
## Run 22 stress 0.2099342 
## ... Procrustes: rmse 0.01542433  max resid 0.152828 
## Run 23 stress 0.2103332 
## Run 24 stress 0.2102011 
## ... Procrustes: rmse 0.006702283  max resid 0.07031149 
## Run 25 stress 0.2202863 
## Run 26 stress 0.2099359 
## ... Procrustes: rmse 0.01602178  max resid 0.152552 
## Run 27 stress 0.2098393 
## ... Procrustes: rmse 0.01134596  max resid 0.1269412 
## Run 28 stress 0.2146908 
## Run 29 stress 0.2169865 
## Run 30 stress 0.2129698 
## Run 31 stress 0.2135891 
## Run 32 stress 0.2170369 
## Run 33 stress 0.2098652 
## ... Procrustes: rmse 0.009361348  max resid 0.1212622 
## Run 34 stress 0.2102694 
## Run 35 stress 0.2159039 
## Run 36 stress 0.2102638 
## Run 37 stress 0.2163209 
## Run 38 stress 0.2103664 
## Run 39 stress 0.2126085 
## Run 40 stress 0.2124092 
## Run 41 stress 0.2129505 
## Run 42 stress 0.2198644 
## Run 43 stress 0.2133866 
## Run 44 stress 0.2254004 
## Run 45 stress 0.2192341 
## Run 46 stress 0.213389 
## Run 47 stress 0.2228373 
## Run 48 stress 0.2099019 
## ... Procrustes: rmse 0.01294875  max resid 0.1525675 
## Run 49 stress 0.209759 
## ... Procrustes: rmse 0.00652526  max resid 0.04372295 
## Run 50 stress 0.2102149 
## ... Procrustes: rmse 0.01700174  max resid 0.1532311 
## Run 51 stress 0.228718 
## Run 52 stress 0.2099745 
## ... Procrustes: rmse 0.01454891  max resid 0.1533141 
## Run 53 stress 0.2130074 
## Run 54 stress 0.2126532 
## Run 55 stress 0.2132173 
## Run 56 stress 0.21023 
## Run 57 stress 0.2136691 
## Run 58 stress 0.2176221 
## Run 59 stress 0.2125307 
## Run 60 stress 0.2130349 
## Run 61 stress 0.2101029 
## ... Procrustes: rmse 0.01033814  max resid 0.09608037 
## Run 62 stress 0.209758 
## ... Procrustes: rmse 0.006722914  max resid 0.04813515 
## Run 63 stress 0.2098131 
## ... Procrustes: rmse 0.008581975  max resid 0.08773394 
## Run 64 stress 0.2278306 
## Run 65 stress 0.2102423 
## Run 66 stress 0.2101894 
## ... Procrustes: rmse 0.00970325  max resid 0.1026161 
## Run 67 stress 0.212608 
## Run 68 stress 0.2168897 
## Run 69 stress 0.2099404 
## ... Procrustes: rmse 0.0163191  max resid 0.1526328 
## Run 70 stress 0.2131248 
## Run 71 stress 0.2122657 
## Run 72 stress 0.2133126 
## Run 73 stress 0.2146229 
## Run 74 stress 0.2099431 
## ... Procrustes: rmse 0.01432381  max resid 0.1532026 
## Run 75 stress 0.2173323 
## Run 76 stress 0.2152484 
## Run 77 stress 0.2117179 
## Run 78 stress 0.2102342 
## Run 79 stress 0.2098498 
## ... Procrustes: rmse 0.008711821  max resid 0.0869956 
## Run 80 stress 0.209855 
## ... Procrustes: rmse 0.0111704  max resid 0.1073192 
## Run 81 stress 0.2130124 
## Run 82 stress 0.2099306 
## ... Procrustes: rmse 0.01604905  max resid 0.1525495 
## Run 83 stress 0.2148915 
## Run 84 stress 0.2133183 
## Run 85 stress 0.2152605 
## Run 86 stress 0.2124853 
## Run 87 stress 0.2112723 
## Run 88 stress 0.2101393 
## ... Procrustes: rmse 0.007157299  max resid 0.07225963 
## Run 89 stress 0.212616 
## Run 90 stress 0.2116955 
## Run 91 stress 0.224748 
## Run 92 stress 0.2097611 
## ... Procrustes: rmse 0.004834871  max resid 0.02894467 
## Run 93 stress 0.2134723 
## Run 94 stress 0.2102688 
## Run 95 stress 0.2098322 
## ... Procrustes: rmse 0.0129102  max resid 0.1366845 
## Run 96 stress 0.2099257 
## ... Procrustes: rmse 0.01588553  max resid 0.1526738 
## Run 97 stress 0.2131946 
## Run 98 stress 0.2160033 
## Run 99 stress 0.2103055 
## Run 100 stress 0.2098845 
## ... Procrustes: rmse 0.0126558  max resid 0.1525857 
## Run 101 stress 0.2098822 
## ... Procrustes: rmse 0.01203187  max resid 0.1530025 
## Run 102 stress 0.209722 
## ... New best solution
## ... Procrustes: rmse 0.005818339  max resid 0.0445689 
## Run 103 stress 0.2133288 
## Run 104 stress 0.2163187 
## Run 105 stress 0.2161437 
## Run 106 stress 0.2123118 
## Run 107 stress 0.2129505 
## Run 108 stress 0.2099305 
## ... Procrustes: rmse 0.01562328  max resid 0.1533819 
## Run 109 stress 0.2097666 
## ... Procrustes: rmse 0.00932821  max resid 0.1214356 
## Run 110 stress 0.2161913 
## Run 111 stress 0.2101285 
## ... Procrustes: rmse 0.007023751  max resid 0.07412396 
## Run 112 stress 0.2097272 
## ... Procrustes: rmse 0.00305532  max resid 0.02186644 
## Run 113 stress 0.2118821 
## Run 114 stress 0.2117277 
## Run 115 stress 0.2098152 
## ... Procrustes: rmse 0.01051393  max resid 0.09773606 
## Run 116 stress 0.2133935 
## Run 117 stress 0.2119973 
## Run 118 stress 0.2129695 
## Run 119 stress 0.2120935 
## Run 120 stress 0.2134367 
## Run 121 stress 0.2126096 
## Run 122 stress 0.2129986 
## Run 123 stress 0.2126035 
## Run 124 stress 0.2097484 
## ... Procrustes: rmse 0.002385626  max resid 0.02715788 
## Run 125 stress 0.2370965 
## Run 126 stress 0.2146537 
## Run 127 stress 0.2099287 
## ... Procrustes: rmse 0.01564831  max resid 0.1533717 
## Run 128 stress 0.2098123 
## ... Procrustes: rmse 0.01282204  max resid 0.1264842 
## Run 129 stress 0.2150643 
## Run 130 stress 0.2099154 
## ... Procrustes: rmse 0.01450937  max resid 0.1534787 
## Run 131 stress 0.2116439 
## Run 132 stress 0.2100327 
## ... Procrustes: rmse 0.01605158  max resid 0.1537851 
## Run 133 stress 0.2115102 
## Run 134 stress 0.2182735 
## Run 135 stress 0.2103294 
## Run 136 stress 0.2098103 
## ... Procrustes: rmse 0.01094055  max resid 0.12632 
## Run 137 stress 0.2166619 
## Run 138 stress 0.2137642 
## Run 139 stress 0.2126141 
## Run 140 stress 0.2127359 
## Run 141 stress 0.210165 
## ... Procrustes: rmse 0.007980591  max resid 0.0748096 
## Run 142 stress 0.2134329 
## Run 143 stress 0.2159204 
## Run 144 stress 0.2101602 
## ... Procrustes: rmse 0.007766581  max resid 0.0746195 
## Run 145 stress 0.2099195 
## ... Procrustes: rmse 0.01187516  max resid 0.153661 
## Run 146 stress 0.2154336 
## Run 147 stress 0.2116886 
## Run 148 stress 0.2101283 
## ... Procrustes: rmse 0.007278672  max resid 0.07457893 
## Run 149 stress 0.2097668 
## ... Procrustes: rmse 0.003413671  max resid 0.03556374 
## Run 150 stress 0.2097016 
## ... New best solution
## ... Procrustes: rmse 0.004209457  max resid 0.02578612 
## Run 151 stress 0.2101298 
## ... Procrustes: rmse 0.0055622  max resid 0.07442458 
## Run 152 stress 0.2112718 
## Run 153 stress 0.2112806 
## Run 154 stress 0.212974 
## Run 155 stress 0.2101334 
## ... Procrustes: rmse 0.007815866  max resid 0.07174455 
## Run 156 stress 0.2097497 
## ... Procrustes: rmse 0.004657534  max resid 0.02896196 
## Run 157 stress 0.2137129 
## Run 158 stress 0.2126415 
## Run 159 stress 0.2097895 
## ... Procrustes: rmse 0.005618096  max resid 0.02922783 
## Run 160 stress 0.2097933 
## ... Procrustes: rmse 0.009126156  max resid 0.1103368 
## Run 161 stress 0.2097696 
## ... Procrustes: rmse 0.006437338  max resid 0.04696808 
## Run 162 stress 0.213439 
## Run 163 stress 0.2103421 
## Run 164 stress 0.2126159 
## Run 165 stress 0.2102967 
## Run 166 stress 0.2121353 
## Run 167 stress 0.2224754 
## Run 168 stress 0.2163261 
## Run 169 stress 0.2248047 
## Run 170 stress 0.210178 
## ... Procrustes: rmse 0.01631992  max resid 0.153225 
## Run 171 stress 0.2097509 
## ... Procrustes: rmse 0.004694005  max resid 0.02844151 
## Run 172 stress 0.2174803 
## Run 173 stress 0.2101205 
## ... Procrustes: rmse 0.005683145  max resid 0.07195479 
## Run 174 stress 0.2150894 
## Run 175 stress 0.2102499 
## Run 176 stress 0.210227 
## Run 177 stress 0.212503 
## Run 178 stress 0.2147534 
## Run 179 stress 0.2326758 
## Run 180 stress 0.2097112 
## ... Procrustes: rmse 0.005422487  max resid 0.04271508 
## Run 181 stress 0.2116783 
## Run 182 stress 0.2098008 
## ... Procrustes: rmse 0.009201616  max resid 0.1118617 
## Run 183 stress 0.2126199 
## Run 184 stress 0.2099307 
## ... Procrustes: rmse 0.01533  max resid 0.1530665 
## Run 185 stress 0.2132181 
## Run 186 stress 0.2098015 
## ... Procrustes: rmse 0.009663798  max resid 0.1227011 
## Run 187 stress 0.2112897 
## Run 188 stress 0.2101601 
## ... Procrustes: rmse 0.01772144  max resid 0.1664074 
## Run 189 stress 0.210365 
## Run 190 stress 0.2097791 
## ... Procrustes: rmse 0.005172674  max resid 0.02794934 
## Run 191 stress 0.2101423 
## ... Procrustes: rmse 0.005995007  max resid 0.07402147 
## Run 192 stress 0.2102844 
## Run 193 stress 0.2098142 
## ... Procrustes: rmse 0.009969578  max resid 0.126371 
## Run 194 stress 0.2130184 
## Run 195 stress 0.2103128 
## Run 196 stress 0.2184558 
## Run 197 stress 0.2102422 
## Run 198 stress 0.213911 
## Run 199 stress 0.2245082 
## Run 200 stress 0.2150762 
## Run 201 stress 0.2225779 
## Run 202 stress 0.2129392 
## Run 203 stress 0.2097176 
## ... Procrustes: rmse 0.005990087  max resid 0.04087518 
## Run 204 stress 0.2132634 
## Run 205 stress 0.2098963 
## ... Procrustes: rmse 0.01236197  max resid 0.1530079 
## Run 206 stress 0.2126056 
## Run 207 stress 0.224634 
## Run 208 stress 0.2147059 
## Run 209 stress 0.2198905 
## Run 210 stress 0.213062 
## Run 211 stress 0.209894 
## ... Procrustes: rmse 0.01226484  max resid 0.1530733 
## Run 212 stress 0.2247818 
## Run 213 stress 0.209932 
## ... Procrustes: rmse 0.01553434  max resid 0.1529955 
## Run 214 stress 0.2136695 
## Run 215 stress 0.2097927 
## ... Procrustes: rmse 0.009248614  max resid 0.1160047 
## Run 216 stress 0.2130952 
## Run 217 stress 0.217344 
## Run 218 stress 0.2125043 
## Run 219 stress 0.2162334 
## Run 220 stress 0.2098528 
## ... Procrustes: rmse 0.01088873  max resid 0.09649557 
## Run 221 stress 0.2099269 
## ... Procrustes: rmse 0.01482054  max resid 0.1528368 
## Run 222 stress 0.2144193 
## Run 223 stress 0.2101275 
## ... Procrustes: rmse 0.005750061  max resid 0.07361562 
## Run 224 stress 0.212064 
## Run 225 stress 0.2146626 
## Run 226 stress 0.2137951 
## Run 227 stress 0.2099296 
## ... Procrustes: rmse 0.013993  max resid 0.1531029 
## Run 228 stress 0.23434 
## Run 229 stress 0.2154847 
## Run 230 stress 0.2228331 
## Run 231 stress 0.2144516 
## Run 232 stress 0.2098424 
## ... Procrustes: rmse 0.01039714  max resid 0.11132 
## Run 233 stress 0.2097181 
## ... Procrustes: rmse 0.00423394  max resid 0.02688357 
## Run 234 stress 0.2145272 
## Run 235 stress 0.2144672 
## Run 236 stress 0.2120379 
## Run 237 stress 0.2161529 
## Run 238 stress 0.2117142 
## Run 239 stress 0.2102623 
## Run 240 stress 0.2161653 
## Run 241 stress 0.2099349 
## ... Procrustes: rmse 0.01458496  max resid 0.1531812 
## Run 242 stress 0.2125044 
## Run 243 stress 0.2126944 
## Run 244 stress 0.2258701 
## Run 245 stress 0.2267006 
## Run 246 stress 0.2134065 
## Run 247 stress 0.2129694 
## Run 248 stress 0.2145001 
## Run 249 stress 0.2123909 
## Run 250 stress 0.2130058 
## Run 251 stress 0.2113364 
## Run 252 stress 0.2099315 
## ... Procrustes: rmse 0.0153496  max resid 0.1530516 
## Run 253 stress 0.2163246 
## Run 254 stress 0.2116937 
## Run 255 stress 0.2099188 
## ... Procrustes: rmse 0.0150663  max resid 0.1530329 
## Run 256 stress 0.2103265 
## Run 257 stress 0.2102338 
## Run 258 stress 0.2097157 
## ... Procrustes: rmse 0.008446499  max resid 0.0961253 
## Run 259 stress 0.2131688 
## Run 260 stress 0.2099365 
## ... Procrustes: rmse 0.01558455  max resid 0.1529922 
## Run 261 stress 0.2101178 
## ... Procrustes: rmse 0.005679899  max resid 0.07294909 
## Run 262 stress 0.2126052 
## Run 263 stress 0.2101146 
## ... Procrustes: rmse 0.005736433  max resid 0.07349935 
## Run 264 stress 0.2097854 
## ... Procrustes: rmse 0.009056932  max resid 0.1018646 
## Run 265 stress 0.2097486 
## ... Procrustes: rmse 0.004576749  max resid 0.02854606 
## Run 266 stress 0.2181454 
## Run 267 stress 0.2101606 
## ... Procrustes: rmse 0.01047877  max resid 0.1044486 
## Run 268 stress 0.2099258 
## ... Procrustes: rmse 0.01433792  max resid 0.1530823 
## Run 269 stress 0.211862 
## Run 270 stress 0.2097949 
## ... Procrustes: rmse 0.006538335  max resid 0.04690468 
## Run 271 stress 0.2102378 
## Run 272 stress 0.2148592 
## Run 273 stress 0.212043 
## Run 274 stress 0.2148504 
## Run 275 stress 0.2116562 
## Run 276 stress 0.2101112 
## ... Procrustes: rmse 0.01043439  max resid 0.1013377 
## Run 277 stress 0.2130211 
## Run 278 stress 0.2103682 
## Run 279 stress 0.2099322 
## ... Procrustes: rmse 0.0154728  max resid 0.1529921 
## Run 280 stress 0.2098838 
## ... Procrustes: rmse 0.01184859  max resid 0.1531203 
## Run 281 stress 0.2148362 
## Run 282 stress 0.2102986 
## Run 283 stress 0.2126003 
## Run 284 stress 0.2129711 
## Run 285 stress 0.2098408 
## ... Procrustes: rmse 0.00992117  max resid 0.122135 
## Run 286 stress 0.2131976 
## Run 287 stress 0.2188947 
## Run 288 stress 0.2125095 
## Run 289 stress 0.2125251 
## Run 290 stress 0.2213913 
## Run 291 stress 0.2103218 
## Run 292 stress 0.2103556 
## Run 293 stress 0.2102009 
## ... Procrustes: rmse 0.006264532  max resid 0.07087182 
## Run 294 stress 0.2128862 
## Run 295 stress 0.2293045 
## Run 296 stress 0.2101116 
## ... Procrustes: rmse 0.005800042  max resid 0.0737227 
## Run 297 stress 0.210243 
## Run 298 stress 0.2134944 
## Run 299 stress 0.2172719 
## Run 300 stress 0.2162612 
## Run 301 stress 0.2434299 
## Run 302 stress 0.2102694 
## Run 303 stress 0.2134215 
## Run 304 stress 0.2127281 
## Run 305 stress 0.2099516 
## ... Procrustes: rmse 0.01221039  max resid 0.1536921 
## Run 306 stress 0.2178865 
## Run 307 stress 0.2097583 
## ... Procrustes: rmse 0.005003836  max resid 0.02866227 
## Run 308 stress 0.2117147 
## Run 309 stress 0.2098077 
## ... Procrustes: rmse 0.009842923  max resid 0.1219967 
## Run 310 stress 0.2165044 
## Run 311 stress 0.212987 
## Run 312 stress 0.2103349 
## Run 313 stress 0.2099146 
## ... Procrustes: rmse 0.01485102  max resid 0.1530597 
## Run 314 stress 0.2169177 
## Run 315 stress 0.2247716 
## Run 316 stress 0.2098862 
## ... Procrustes: rmse 0.01220158  max resid 0.1531984 
## Run 317 stress 0.2133246 
## Run 318 stress 0.2167595 
## Run 319 stress 0.2098052 
## ... Procrustes: rmse 0.00953829  max resid 0.1244171 
## Run 320 stress 0.2112655 
## Run 321 stress 0.2116529 
## Run 322 stress 0.2225972 
## Run 323 stress 0.2102479 
## Run 324 stress 0.2113148 
## Run 325 stress 0.2175125 
## Run 326 stress 0.2120987 
## Run 327 stress 0.210256 
## Run 328 stress 0.2134304 
## Run 329 stress 0.2102517 
## Run 330 stress 0.2129867 
## Run 331 stress 0.2150343 
## Run 332 stress 0.2098352 
## ... Procrustes: rmse 0.00919755  max resid 0.09367652 
## Run 333 stress 0.2168579 
## Run 334 stress 0.2118603 
## Run 335 stress 0.2101493 
## ... Procrustes: rmse 0.01068699  max resid 0.1049335 
## Run 336 stress 0.2102338 
## Run 337 stress 0.2099255 
## ... Procrustes: rmse 0.01505611  max resid 0.1528571 
## Run 338 stress 0.2099145 
## ... Procrustes: rmse 0.0116063  max resid 0.0973198 
## Run 339 stress 0.2099309 
## ... Procrustes: rmse 0.01533235  max resid 0.1530627 
## Run 340 stress 0.2097275 
## ... Procrustes: rmse 0.00395705  max resid 0.02800491 
## Run 341 stress 0.2153898 
## Run 342 stress 0.2394691 
## Run 343 stress 0.21273 
## Run 344 stress 0.2102831 
## Run 345 stress 0.2097241 
## ... Procrustes: rmse 0.004416164  max resid 0.0280246 
## Run 346 stress 0.2098928 
## ... Procrustes: rmse 0.01240304  max resid 0.152975 
## Run 347 stress 0.2097859 
## ... Procrustes: rmse 0.005706258  max resid 0.06201415 
## Run 348 stress 0.2134563 
## Run 349 stress 0.2121168 
## Run 350 stress 0.210378 
## Run 351 stress 0.2103764 
## Run 352 stress 0.2130106 
## Run 353 stress 0.2170193 
## Run 354 stress 0.2098083 
## ... Procrustes: rmse 0.005275467  max resid 0.02941205 
## Run 355 stress 0.2119085 
## Run 356 stress 0.2103907 
## Run 357 stress 0.2102628 
## Run 358 stress 0.2101609 
## ... Procrustes: rmse 0.005813887  max resid 0.07415652 
## Run 359 stress 0.2099309 
## ... Procrustes: rmse 0.01456246  max resid 0.1532625 
## Run 360 stress 0.2175221 
## Run 361 stress 0.2217646 
## Run 362 stress 0.2102498 
## Run 363 stress 0.2102884 
## Run 364 stress 0.2169011 
## Run 365 stress 0.2132479 
## Run 366 stress 0.2099212 
## ... Procrustes: rmse 0.01513998  max resid 0.1529961 
## Run 367 stress 0.2131867 
## Run 368 stress 0.2131265 
## Run 369 stress 0.2116568 
## Run 370 stress 0.2231066 
## Run 371 stress 0.2098369 
## ... Procrustes: rmse 0.009911654  max resid 0.1142778 
## Run 372 stress 0.2097854 
## ... Procrustes: rmse 0.008789799  max resid 0.1177037 
## Run 373 stress 0.2097412 
## ... Procrustes: rmse 0.007941668  max resid 0.08885918 
## Run 374 stress 0.2102052 
## Run 375 stress 0.2097891 
## ... Procrustes: rmse 0.006550043  max resid 0.05198666 
## Run 376 stress 0.213326 
## Run 377 stress 0.2097025 
## ... Procrustes: rmse 0.002876355  max resid 0.01595654 
## Run 378 stress 0.2099077 
## ... Procrustes: rmse 0.01204289  max resid 0.1534777 
## Run 379 stress 0.212503 
## Run 380 stress 0.2148804 
## Run 381 stress 0.2131164 
## Run 382 stress 0.2141244 
## Run 383 stress 0.2098294 
## ... Procrustes: rmse 0.01036347  max resid 0.1322424 
## Run 384 stress 0.2102495 
## Run 385 stress 0.2102679 
## Run 386 stress 0.2179382 
## Run 387 stress 0.2198484 
## Run 388 stress 0.2097184 
## ... Procrustes: rmse 0.004210321  max resid 0.02849628 
## Run 389 stress 0.212339 
## Run 390 stress 0.2148237 
## Run 391 stress 0.2099264 
## ... Procrustes: rmse 0.01534003  max resid 0.1529914 
## Run 392 stress 0.2102592 
## Run 393 stress 0.2101307 
## ... Procrustes: rmse 0.005484187  max resid 0.07336651 
## Run 394 stress 0.209805 
## ... Procrustes: rmse 0.009558297  max resid 0.120224 
## Run 395 stress 0.2163558 
## Run 396 stress 0.2102291 
## Run 397 stress 0.2161451 
## Run 398 stress 0.2126068 
## Run 399 stress 0.2099261 
## ... Procrustes: rmse 0.01509752  max resid 0.1528582 
## Run 400 stress 0.2347234 
## Run 401 stress 0.2112807 
## Run 402 stress 0.2097053 
## ... Procrustes: rmse 0.004117119  max resid 0.02287897 
## Run 403 stress 0.2107541 
## Run 404 stress 0.209775 
## ... Procrustes: rmse 0.009373035  max resid 0.1202368 
## Run 405 stress 0.2151085 
## Run 406 stress 0.237914 
## Run 407 stress 0.2131815 
## Run 408 stress 0.213702 
## Run 409 stress 0.2277012 
## Run 410 stress 0.2097137 
## ... Procrustes: rmse 0.003256294  max resid 0.02391741 
## Run 411 stress 0.2126363 
## Run 412 stress 0.2098564 
## ... Procrustes: rmse 0.009359887  max resid 0.09464593 
## Run 413 stress 0.2096806 
## ... New best solution
## ... Procrustes: rmse 0.003771369  max resid 0.02563352 
## Run 414 stress 0.2101414 
## ... Procrustes: rmse 0.01169878  max resid 0.1053441 
## Run 415 stress 0.211258 
## Run 416 stress 0.21295 
## Run 417 stress 0.2146358 
## Run 418 stress 0.2099262 
## ... Procrustes: rmse 0.01465429  max resid 0.153668 
## Run 419 stress 0.2097208 
## ... Procrustes: rmse 0.009068168  max resid 0.1005321 
## Run 420 stress 0.2099474 
## ... Procrustes: rmse 0.01508528  max resid 0.1535732 
## Run 421 stress 0.2097157 
## ... Procrustes: rmse 0.004738058  max resid 0.04168663 
## Run 422 stress 0.2117891 
## Run 423 stress 0.2179079 
## Run 424 stress 0.2099126 
## ... Procrustes: rmse 0.01355478  max resid 0.1535678 
## Run 425 stress 0.2102493 
## Run 426 stress 0.2220886 
## Run 427 stress 0.2099387 
## ... Procrustes: rmse 0.01348499  max resid 0.1537097 
## Run 428 stress 0.2098314 
## ... Procrustes: rmse 0.0076032  max resid 0.09563233 
## Run 429 stress 0.2112498 
## Run 430 stress 0.2101959 
## Run 431 stress 0.2134495 
## Run 432 stress 0.2134059 
## Run 433 stress 0.2098668 
## ... Procrustes: rmse 0.0103811  max resid 0.08894867 
## Run 434 stress 0.2226698 
## Run 435 stress 0.2130046 
## Run 436 stress 0.2098851 
## ... Procrustes: rmse 0.01242216  max resid 0.1534893 
## Run 437 stress 0.2157712 
## Run 438 stress 0.2166472 
## Run 439 stress 0.2180417 
## Run 440 stress 0.2097523 
## ... Procrustes: rmse 0.003115454  max resid 0.02803743 
## Run 441 stress 0.2148538 
## Run 442 stress 0.2103 
## Run 443 stress 0.2117119 
## Run 444 stress 0.2097551 
## ... Procrustes: rmse 0.006921319  max resid 0.08624046 
## Run 445 stress 0.2099123 
## ... Procrustes: rmse 0.0138407  max resid 0.1535531 
## Run 446 stress 0.2097192 
## ... Procrustes: rmse 0.005393117  max resid 0.04242808 
## Run 447 stress 0.2148501 
## Run 448 stress 0.2195586 
## Run 449 stress 0.2096752 
## ... New best solution
## ... Procrustes: rmse 0.00120412  max resid 0.01553874 
## Run 450 stress 0.2099407 
## ... Procrustes: rmse 0.01512919  max resid 0.1536836 
## Run 451 stress 0.2177598 
## Run 452 stress 0.2133447 
## Run 453 stress 0.2097955 
## ... Procrustes: rmse 0.003885868  max resid 0.02744706 
## Run 454 stress 0.2097694 
## ... Procrustes: rmse 0.005719434  max resid 0.04377941 
## Run 455 stress 0.2101171 
## ... Procrustes: rmse 0.01020592  max resid 0.09090328 
## Run 456 stress 0.209809 
## ... Procrustes: rmse 0.004514653  max resid 0.02835068 
## Run 457 stress 0.2102806 
## Run 458 stress 0.2099179 
## ... Procrustes: rmse 0.01418331  max resid 0.1533645 
## Run 459 stress 0.2126038 
## Run 460 stress 0.2101389 
## ... Procrustes: rmse 0.008306194  max resid 0.07348261 
## Run 461 stress 0.2220578 
## Run 462 stress 0.2096898 
## ... Procrustes: rmse 0.002116752  max resid 0.010882 
## Run 463 stress 0.2120487 
## Run 464 stress 0.2103361 
## Run 465 stress 0.2194496 
## Run 466 stress 0.2175507 
## Run 467 stress 0.2098872 
## ... Procrustes: rmse 0.0119993  max resid 0.1534727 
## Run 468 stress 0.2132519 
## Run 469 stress 0.2129907 
## Run 470 stress 0.2102533 
## Run 471 stress 0.212572 
## Run 472 stress 0.2102243 
## Run 473 stress 0.2112709 
## Run 474 stress 0.2349923 
## Run 475 stress 0.2126826 
## Run 476 stress 0.2154687 
## Run 477 stress 0.2112593 
## Run 478 stress 0.2102025 
## Run 479 stress 0.2147848 
## Run 480 stress 0.2097926 
## ... Procrustes: rmse 0.005780708  max resid 0.04216166 
## Run 481 stress 0.2146509 
## Run 482 stress 0.2260598 
## Run 483 stress 0.2116651 
## Run 484 stress 0.2098007 
## ... Procrustes: rmse 0.004470413  max resid 0.02955439 
## Run 485 stress 0.2099029 
## ... Procrustes: rmse 0.01242719  max resid 0.1534549 
## Run 486 stress 0.2102856 
## Run 487 stress 0.2116572 
## Run 488 stress 0.2124048 
## Run 489 stress 0.211685 
## Run 490 stress 0.2096751 
## ... New best solution
## ... Procrustes: rmse 0.0005534366  max resid 0.004358859 
## ... Similar to previous best
## *** Solution reached
str(meta.nmds.FloristicSurvey2D)
## List of 35
##  $ nobj      : int 196
##  $ nfix      : int 0
##  $ ndim      : num 2
##  $ ndis      : int 19110
##  $ ngrp      : int 1
##  $ diss      : num [1:19110] 0 0 0 0 0 0 0 0 0 0 ...
##  $ iidx      : int [1:19110] 117 111 67 129 186 123 126 129 129 170 ...
##  $ jidx      : int [1:19110] 106 103 2 122 128 121 122 127 126 3 ...
##  $ xinit     : num [1:392] 0.849 0.803 0.633 0.357 0.187 ...
##  $ istart    : int 1
##  $ isform    : int 1
##  $ ities     : int 1
##  $ iregn     : int 1
##  $ iscal     : int 1
##  $ maxits    : int 200
##  $ sratmx    : num 1
##  $ strmin    : num 1e-04
##  $ sfgrmn    : num 1e-07
##  $ dist      : num [1:19110] 0 0 0 0 0 0 0 0 0 0 ...
##  $ dhat      : num [1:19110] 0 0 0 0 0 0 0 0 0 0 ...
##  $ points    : num [1:196, 1:2] -0.136 -0.149 -0.154 -0.105 -0.103 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:196] "7" "7" "7" "7" ...
##   .. ..$ : chr [1:2] "MDS1" "MDS2"
##   ..- attr(*, "centre")= logi TRUE
##   ..- attr(*, "pc")= logi TRUE
##   ..- attr(*, "halfchange")= logi FALSE
##   ..- attr(*, "internalscaling")= num 4.31
##  $ stress    : num 0.21
##  $ grstress  : num 0.21
##  $ iters     : int 148
##  $ icause    : int 3
##  $ call      : language metaMDS(comm = dist_FloristicSurvey, k = 2, trymax = 1000)
##  $ model     : chr "global"
##  $ distmethod: chr "binary bray"
##  $ distcall  : chr "vegdist(x = M, method = \"bray\", binary = T)"
##  $ distance  : chr "binary bray"
##  $ converged : logi TRUE
##  $ tries     : num 490
##  $ engine    : chr "monoMDS"
##  $ species   : logi NA
##  $ data      : chr "dist_FloristicSurvey"
##  - attr(*, "class")= chr [1:2] "metaMDS" "monoMDS"
stressplot(meta.nmds.FloristicSurvey2D)

FloristicSurvey_envfit <- envfit(meta.nmds.FloristicSurvey2D, env = GM_coverage_df, perm = 999) #standard envfit
FloristicSurvey_envfit
## 
## ***FACTORS:
## 
## Centroids:
##                   NMDS1   NMDS2
## M...52.High      0.0453  0.0401
## M...52.Low      -0.0896 -0.0251
## M...52.Med_High  0.0782 -0.0018
## M...52.Med_Low   0.0455  0.0313
## 
## Goodness of fit:
##             r2 Pr(>r)    
## M...52. 0.1181  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
#data for plotting 
##NMDS points
FloristicSurvey.NMDS.data<-GM_coverage_df 
FloristicSurvey.NMDS.data$NMDS1<-meta.nmds.FloristicSurvey2D$points[ ,1] 
FloristicSurvey.NMDS.data$NMDS2<-meta.nmds.FloristicSurvey2D$points[ ,2] 
colnames(FloristicSurvey.NMDS.data)[1] <- "GM_Coverage"

# data for the envfit arrows
env.scores.FloristicSurvey <- as.data.frame(scores(FloristicSurvey_envfit, display = "vectors")) #extracts relevant scores from envifit
env.scores.FloristicSurvey <- cbind(env.scores.FloristicSurvey, env.variables = rownames(env.scores.FloristicSurvey)) #and then gives them their names

# function for ellipsess - just run this, is used later
#taken from the excellent stackoverflow Q+A: http://stackoverflow.com/questions/13794419/plotting-ordiellipse-function-from-vegan-package-onto-nmds-plot-created-in-ggplo
veganCovEllipse <- function (cov, center = c(0, 0), scale = 1, npoints = 100) 
{
  theta <- (0:npoints) * 2 * pi/npoints
  Circle <- cbind(cos(theta), sin(theta))
  t(center + scale * t(Circle %*% chol(cov)))
}

#data for ellipse, use GM coverage
df_ell.FSurvey.GM_coverage <- data.frame() #sets up a data frame before running the function.
for(g in levels(FloristicSurvey.NMDS.data$GM_Coverage)){
  df_ell.FSurvey.GM_coverage <- rbind(df_ell.FSurvey.GM_coverage, cbind(as.data.frame(with(FloristicSurvey.NMDS.data [FloristicSurvey.NMDS.data$GM_Coverage==g,],
                                                                                   veganCovEllipse(cov.wt(cbind(NMDS1,NMDS2),wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS2)))))
                                                                ,GM_coverage=g))
}

# data for labelling the ellipse
NMDS.mean.FloristicSurvey=aggregate(FloristicSurvey.NMDS.data[ ,c("NMDS1", "NMDS2")], 
                         list(group = FloristicSurvey.NMDS.data$GM_Coverage), mean)

## finally plotting.
mult <- 1 #multiplier for the arrows and text for envfit below. You can change this and then rerun the plot command.
ggplot(data = FloristicSurvey.NMDS.data, aes(y = NMDS2, x = NMDS1))+ 
    geom_path(data = df_ell.FSurvey.GM_coverage, aes(x = NMDS1, y = NMDS2, group = df_ell.FSurvey.GM_coverage$GM_coverage, color=df_ell.FSurvey.GM_coverage$GM_coverage))+ 
    geom_point( aes(color = FloristicSurvey.NMDS.data$GM_Coverage), size = 1)+ 
    scale_color_manual(values = c(14,5,19,2))+ 
    coord_cartesian(xlim = c(-1,1.5))+
  theme_simple()+ guides(color=guide_legend("Garlic Mustard Coverage"))

PCA

FloristicSurvey<-FloristicSurvey[,-56]
FloristicSurvey.pca <- prcomp(FloristicSurvey[,4:(ncol(FloristicSurvey))],center = TRUE) 
summary(FloristicSurvey.pca)
## Importance of components%s:
##                           PC1    PC2     PC3     PC4     PC5    PC6
## Standard deviation     0.6406 0.5639 0.49314 0.45065 0.42882 0.3954
## Proportion of Variance 0.1365 0.1058 0.08088 0.06755 0.06116 0.0520
## Cumulative Proportion  0.1365 0.2422 0.32311 0.39066 0.45182 0.5038
##                            PC7    PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.37368 0.3608 0.34564 0.33141 0.30059 0.27972
## Proportion of Variance 0.04644 0.0433 0.03974 0.03653 0.03005 0.02602
## Cumulative Proportion  0.55027 0.5936 0.63330 0.66983 0.69989 0.72591
##                           PC13    PC14    PC15    PC16    PC17   PC18
## Standard deviation     0.26731 0.23707 0.23511 0.23170 0.23140 0.2159
## Proportion of Variance 0.02377 0.01869 0.01839 0.01786 0.01781 0.0155
## Cumulative Proportion  0.74968 0.76837 0.78675 0.80461 0.82242 0.8379
##                           PC19   PC20   PC21    PC22    PC23    PC24
## Standard deviation     0.20362 0.1992 0.1835 0.17727 0.17130 0.16138
## Proportion of Variance 0.01379 0.0132 0.0112 0.01045 0.00976 0.00866
## Cumulative Proportion  0.85171 0.8649 0.8761 0.88656 0.89632 0.90499
##                           PC25    PC26    PC27    PC28    PC29    PC30
## Standard deviation     0.16000 0.15304 0.14850 0.14015 0.13172 0.12966
## Proportion of Variance 0.00851 0.00779 0.00733 0.00653 0.00577 0.00559
## Cumulative Proportion  0.91350 0.92129 0.92863 0.93516 0.94093 0.94652
##                           PC31    PC32    PC33    PC34    PC35    PC36
## Standard deviation     0.12335 0.11774 0.11421 0.11329 0.10677 0.10455
## Proportion of Variance 0.00506 0.00461 0.00434 0.00427 0.00379 0.00364
## Cumulative Proportion  0.95158 0.95619 0.96053 0.96480 0.96859 0.97223
##                           PC37    PC38    PC39    PC40    PC41    PC42
## Standard deviation     0.10160 0.09716 0.09407 0.08205 0.07828 0.07736
## Proportion of Variance 0.00343 0.00314 0.00294 0.00224 0.00204 0.00199
## Cumulative Proportion  0.97566 0.97880 0.98174 0.98398 0.98602 0.98801
##                           PC43    PC44    PC45    PC46    PC47    PC48
## Standard deviation     0.07252 0.06984 0.06822 0.06603 0.06268 0.05955
## Proportion of Variance 0.00175 0.00162 0.00155 0.00145 0.00131 0.00118
## Cumulative Proportion  0.98976 0.99138 0.99293 0.99438 0.99569 0.99687
##                           PC49    PC50    PC51      PC52
## Standard deviation     0.05928 0.05691 0.05170 4.972e-17
## Proportion of Variance 0.00117 0.00108 0.00089 0.000e+00
## Cumulative Proportion  0.99803 0.99911 1.00000 1.000e+00
fviz_eig(FloristicSurvey.pca,addlabels=TRUE,hjust = -0.3)+theme_simple()

round(FloristicSurvey.pca$sdev^2/sum(FloristicSurvey.pca$sdev^2)*100,2)
##  [1] 13.65 10.58  8.09  6.75  6.12  5.20  4.64  4.33  3.97  3.65  3.01
## [12]  2.60  2.38  1.87  1.84  1.79  1.78  1.55  1.38  1.32  1.12  1.05
## [23]  0.98  0.87  0.85  0.78  0.73  0.65  0.58  0.56  0.51  0.46  0.43
## [34]  0.43  0.38  0.36  0.34  0.31  0.29  0.22  0.20  0.20  0.17  0.16
## [45]  0.15  0.15  0.13  0.12  0.12  0.11  0.09  0.00
# About 25% explained by my first 2 eigenvectors. 50% by the first 6.

Modeling (redundancy analysis)

species comp ~ GMcov * pop + soil characteristics inside and outside as covariates?

species<-names(FloristicSurvey[4:(length(FloristicSurvey))])
SpeciesGMCorrelation <- data.frame(matrix(ncol = 3, nrow = 0))
names(SpeciesGMCorrelation)<-c("Species", "PValue","signif")

for(i in 4:(length(FloristicSurvey))){
  SpeciesGMCorrelation[i-3,1]<-names(FloristicSurvey[i])
  SpeciesGMCorrelation[i-3,2]<-summary(glm(FloristicSurvey[,i]~GM_Coverage, data=FloristicSurvey, family = binomial))$coefficients[8]
}
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
for(i in 1:length(species)) if( SpeciesGMCorrelation[i,2] <= 0.05){
  SpeciesGMCorrelation[i,3]<-c("TRUE")
}else{
  SpeciesGMCorrelation[i,3]<-c("FALSE")
  }
length(SpeciesGMCorrelation$signif[SpeciesGMCorrelation[,3]==TRUE])/length(SpeciesGMCorrelation$signif)*100
## [1] 25
SpeciesGMCorrelation$Species[SpeciesGMCorrelation[,3]==TRUE]
##  [1] "Acer_Saccharum"              "Galium_circaezans"          
##  [3] "Ulmus_americana"             "Aliaria_petiolata"          
##  [5] "Parthenocissus_quinquefolia" "Rhubarb"                    
##  [7] "Trifolium_sp"                "Rubus_occidentalis"         
##  [9] "Circeae_lutetiana"           "Clematis_virginiana"        
## [11] "Rhamnus_cathartica"          "Unknown_8"                  
## [13] "unknown_12"
# Redundancy analysis
sp.comp <- FloristicSurvey[,4:(length(FloristicSurvey))]
var <- FloristicSurvey[,1:3]
var$Quadrat <- as.character(var$Quadrat)
for(i in 1:length(var$Quadrat))if(var$Population[i] == "7"){
  var$Quadrat[i]<-paste("7",var$Quadrat[i],sep="")
}else if(var$Population[i] == "3"){
  var$Quadrat[i]<-paste("3",var$Quadrat[i],sep="")
}else if(var$Population[i] == "14"){
  var$Quadrat[i]<-paste("14",var$Quadrat[i],sep="")
}else if(var$Population[i] == "13"){
  var$Quadrat[i]<-paste("13",var$Quadrat[i],sep="")
}else if(var$Population[i] == "1"){
  var$Quadrat[i]<-paste("1",var$Quadrat[i],sep="")
}

var$Population <- as.factor(var$Population)

formulaRDA <- rda(sp.comp ~ GM_Coverage + Condition(Population) , data=var)

formulaRDA <- rda(sp.comp ~ GM_Coverage , data=var)
head(summary(formulaRDA))
## 
## Call:
## rda(formula = sp.comp ~ GM_Coverage, data = var) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          3.0066    1.00000
## Constrained    0.1451    0.04825
## Unconstrained  2.8615    0.95175
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          RDA1    PC1     PC2     PC3     PC4     PC5
## Eigenvalue            0.14507 0.3830 0.25896 0.23497 0.19681 0.18378
## Proportion Explained  0.04825 0.1274 0.08613 0.07815 0.06546 0.06113
## Cumulative Proportion 0.04825 0.1757 0.26178 0.33993 0.40539 0.46652
##                           PC6    PC7     PC8     PC9    PC10    PC11
## Eigenvalue            0.14105 0.1302 0.12555 0.11261 0.10979 0.09025
## Proportion Explained  0.04691 0.0433 0.04176 0.03746 0.03652 0.03002
## Cumulative Proportion 0.51343 0.5567 0.59849 0.63595 0.67246 0.70248
##                          PC12    PC13    PC14    PC15    PC16    PC17
## Eigenvalue            0.07818 0.07134 0.05604 0.05392 0.05355 0.05152
## Proportion Explained  0.02600 0.02373 0.01864 0.01793 0.01781 0.01714
## Cumulative Proportion 0.72848 0.75221 0.77085 0.78878 0.80659 0.82373
##                          PC18    PC19   PC20    PC21    PC22    PC23
## Eigenvalue            0.04588 0.04146 0.0394 0.03355 0.03051 0.02929
## Proportion Explained  0.01526 0.01379 0.0131 0.01116 0.01015 0.00974
## Cumulative Proportion 0.83899 0.85278 0.8659 0.87704 0.88719 0.89693
##                          PC24    PC25    PC26    PC27    PC28    PC29
## Eigenvalue            0.02567 0.02559 0.02324 0.02183 0.01964 0.01731
## Proportion Explained  0.00854 0.00851 0.00773 0.00726 0.00653 0.00576
## Cumulative Proportion 0.90547 0.91398 0.92171 0.92897 0.93551 0.94126
##                          PC30    PC31    PC32    PC33    PC34    PC35
## Eigenvalue            0.01681 0.01508 0.01386 0.01304 0.01283 0.01139
## Proportion Explained  0.00559 0.00501 0.00461 0.00434 0.00427 0.00379
## Cumulative Proportion 0.94686 0.95187 0.95648 0.96082 0.96508 0.96887
##                          PC36    PC37     PC38     PC39     PC40     PC41
## Eigenvalue            0.01089 0.01027 0.009383 0.008708 0.006676 0.006114
## Proportion Explained  0.00362 0.00342 0.003120 0.002900 0.002220 0.002030
## Cumulative Proportion 0.97249 0.97591 0.979030 0.981930 0.984150 0.986180
##                           PC42     PC43     PC44     PC45     PC46
## Eigenvalue            0.005879 0.005254 0.004851 0.004653 0.004186
## Proportion Explained  0.001960 0.001750 0.001610 0.001550 0.001390
## Cumulative Proportion 0.988140 0.989880 0.991500 0.993050 0.994440
##                           PC47     PC48     PC49     PC50     PC51
## Eigenvalue            0.003919 0.003546 0.003458 0.003175 0.002623
## Proportion Explained  0.001300 0.001180 0.001150 0.001060 0.000870
## Cumulative Proportion 0.995740 0.996920 0.998070 0.999130 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1
## Eigenvalue            0.1451
## Proportion Explained  1.0000
## Cumulative Proportion 1.0000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  4.933262 
## 
## 
## Species scores
## 
##                        RDA1      PC1      PC2       PC3      PC4      PC5
## Galium_aparine     -0.03884  0.95611 -0.68111  0.398026  0.56161 -0.16776
## Acer_Saccharum     -0.39518 -0.59225 -0.76442 -0.724167  0.39649  0.40523
## Unknown_Grass      -0.09710  0.83165  0.29277 -0.886732 -0.01378 -0.05976
## Galium_circaezans  -0.09565  0.08758 -0.06161 -0.009061  0.03333 -0.06610
## Ulmus_americana    -0.18939  0.19041 -0.59667 -0.292853 -0.51085 -0.72307
## False_Solomon_Seal -0.05895  0.05660  0.00144 -0.025574  0.01268 -0.01885
## ....                                                                     
## 
## 
## Site scores (weighted sums of species scores)
## 
##         RDA1    PC1     PC2     PC3       PC4      PC5
## row1 -0.8792 0.3117 -0.4541 -0.3298  0.088931 -0.59497
## row2 -0.9951 0.2201 -0.4529 -0.3803 -0.087225 -0.39023
## row3 -0.6743 0.1376  0.2111 -0.5348 -0.312633  0.51839
## row4 -0.8687 0.5073 -0.3985 -0.4795 -0.280767 -0.08044
## row5 -0.7576 0.4642 -0.1132 -0.3390  0.047862  0.39295
## row6 -0.5475 0.3770  0.2156  0.1561  0.003758 -0.18769
## ....                                                  
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         RDA1    PC1     PC2     PC3       PC4      PC5
## row1 -0.4053 0.3117 -0.4541 -0.3298  0.088931 -0.59497
## row2 -0.4053 0.2201 -0.4529 -0.3803 -0.087225 -0.39023
## row3 -0.4591 0.1376  0.2111 -0.5348 -0.312633  0.51839
## row4 -0.4591 0.5073 -0.3985 -0.4795 -0.280767 -0.08044
## row5 -0.4591 0.4642 -0.1132 -0.3390  0.047862  0.39295
## row6 -0.4591 0.3770  0.2156  0.1561  0.003758 -0.18769
## ....                                                  
## 
## 
## Biplot scores for constraining variables
## 
##             RDA1 PC1 PC2 PC3 PC4 PC5
## GM_Coverage    1   0   0   0   0   0
RsquareAdj(formulaRDA)
## $r.squared
## [1] 0.04825094
## 
## $adj.r.squared
## [1] 0.04339508
plot(formulaRDA, scaling = 2)

anova.cca(formulaRDA,step=1000)
anova.cca(formulaRDA,by="axis", step=1000)

*Note that I redid this with GM_coverage_catogory and i didnt find ANY significance!

Soil Characteristics

Soil_characteristics$Location<-as.factor(Soil_characteristics$Location)
GLmodel_soil <- glm(Location ~Agg_stability*C_percent*N_percent,family=binomial(link='logit'),data=Soil_characteristics)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(GLmodel_soil)
## 
## Call:
## glm(formula = Location ~ Agg_stability * C_percent * N_percent, 
##     family = binomial(link = "logit"), data = Soil_characteristics)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5019  -0.9103  -0.1493   0.7927   2.1777  
## 
## Coefficients:
##                                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)                         -40.68      40.25  -1.011    0.312
## Agg_stability                        57.36      62.43   0.919    0.358
## C_percent                            28.40      18.23   1.558    0.119
## N_percent                          -178.24     138.30  -1.289    0.197
## Agg_stability:C_percent             -40.73      26.92  -1.513    0.130
## Agg_stability:N_percent             253.80     211.09   1.202    0.229
## C_percent:N_percent                 -14.38      10.04  -1.432    0.152
## Agg_stability:C_percent:N_percent    20.95      15.18   1.380    0.168
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 27.726  on 19  degrees of freedom
## Residual deviance: 18.370  on 12  degrees of freedom
## AIC: 34.37
## 
## Number of Fisher Scoring iterations: 8
Soil_characteristics2 <- data.frame(Soil_characteristics$Population, Soil_characteristics$Location, log(Soil_characteristics$pH), log(Soil_characteristics$Agg_stability), log(Soil_characteristics$C_percent), log(Soil_characteristics$N_percent))
Soil_characteristics2<-rename(Soil_characteristics2, c("Soil_characteristics.Population"="Population", "Soil_characteristics.Location"="Location", "log.Soil_characteristics.pH."="pH", "log.Soil_characteristics.Agg_stability."="Agg_stability", "log.Soil_characteristics.C_percent."="Percent","log.Soil_characteristics.N_percent."="Percent" ))
## Warning: The plyr::rename operation has created duplicates for the
## following name(s): (`Percent`)
df2 <- reshape2::melt(Soil_characteristics2, id.var=c("Population","Location"))
df2 <- rename(df2, c("value"="Log_value"))

par(mfrow=c(2,2))
ggplot(Soil_characteristics, aes(Location,Agg_stability,fill=Location))+ geom_boxplot() +
  theme_simple()

ggplot(Soil_characteristics, aes(Location,pH,fill=Location))+ geom_boxplot() +
  theme_simple()

ggplot(Soil_characteristics, aes(Location,C_percent,fill=Location))+ geom_boxplot() +
  theme_simple()

ggplot(Soil_characteristics, aes(Location,N_percent,fill=Location))+ geom_boxplot() +
  theme_simple()

Soil_characteristics$Population<-as.factor(Soil_characteristics$Population)

GLMER_soil <- glmer(Location~scale(Agg_stability)*scale(C_percent)*scale(N_percent)+(1|Population),family=binomial(link='logit'),data=Soil_characteristics)
summary(GLMER_soil)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Location ~ scale(Agg_stability) * scale(C_percent) * scale(N_percent) +  
##     (1 | Population)
##    Data: Soil_characteristics
## 
##      AIC      BIC   logLik deviance df.resid 
##     36.4     45.3     -9.2     18.4       11 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4454 -0.7165 -0.1067  0.6076  3.1160 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev. 
##  Population (Intercept) 6.629e-17 8.142e-09
## Number of obs: 20, groups:  Population, 10
## 
## Fixed effects:
##                                                        Estimate Std. Error
## (Intercept)                                               1.895      1.494
## scale(Agg_stability)                                     -3.634      2.731
## scale(C_percent)                                         10.664      7.517
## scale(N_percent)                                        -10.010      7.247
## scale(Agg_stability):scale(C_percent)                   -15.686     10.695
## scale(Agg_stability):scale(N_percent)                    15.667     10.403
## scale(C_percent):scale(N_percent)                        -1.937      1.343
## scale(Agg_stability):scale(C_percent):scale(N_percent)    3.464      2.511
##                                                        z value Pr(>|z|)
## (Intercept)                                              1.269    0.204
## scale(Agg_stability)                                    -1.331    0.183
## scale(C_percent)                                         1.419    0.156
## scale(N_percent)                                        -1.381    0.167
## scale(Agg_stability):scale(C_percent)                   -1.467    0.142
## scale(Agg_stability):scale(N_percent)                    1.506    0.132
## scale(C_percent):scale(N_percent)                       -1.442    0.149
## scale(Agg_stability):scale(C_percent):scale(N_percent)   1.380    0.168
## 
## Correlation of Fixed Effects:
##             (Intr) sc(A_) sc(C_) sc(N_) sc(A_):(C_) s(A_):(N s(C_):
## scl(Agg_st) -0.740                                                 
## scl(C_prcn)  0.617 -0.441                                          
## scl(N_prcn) -0.522  0.384 -0.986                                   
## sc(A_):(C_) -0.608  0.791 -0.637  0.604                            
## sc(A_):(N_)  0.558 -0.700  0.635 -0.602 -0.981                     
## sc(C_):(N_) -0.829  0.730 -0.454  0.350  0.695      -0.678         
## s(A_):(C_):  0.718 -0.866  0.656 -0.629 -0.774       0.684   -0.635
options(na.action = "na.fail")
dredge(GLMER_soil) 
## Fixed term is "(Intercept)"
# A model with nothing is best...

Histology data Analysis

Histology Data exploration

Calculate plant-level Mycchorizae, pathogens and herbivory:

# Fix column characteristics for plotting and analysis
fdata$indiv2<-as.factor(fdata$indiv2)
fdata$Mycorrhiza<-as.numeric(paste(fdata$Mycorrhiza))
fdata$Lesion<-as.numeric(paste(fdata$Lesion))
fdata$Herbivory<-as.numeric(paste(fdata$Herbivory))
pct<-function(x){
  sum(x)/100
}

# Calculate plant-level averages
fdat_ind<-aggregate(fdata[,c("Mycorrhiza","Lesion","Herbivory")],by=list(indiv2=fdata$indiv2,pop=fdata$pop,loc=fdata$loc,species=fdata$species), FUN=pct)

Simple linear model visualization:

## Mycorrhiza model (NOTE: No loc*species because not all species present in all populations)
base.Myc<-lmer(Mycorrhiza~species+loc+species:loc+(1|pop)+(1|pop:loc),data=fdat_ind)
noint.Myc<-lmer(Mycorrhiza~species+loc+(1|pop)+(1|pop:loc),data=fdat_ind)
anova(base.Myc,noint.Myc)
## refitting model(s) with ML (instead of REML)
base.Myc ## Inspect coefficients
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## Mycorrhiza ~ species + loc + species:loc + (1 | pop) + (1 | pop:loc)
##    Data: fdat_ind
## REML criterion at convergence: -12.5748
## Random effects:
##  Groups   Name        Std.Dev.
##  pop:loc  (Intercept) 0.06622 
##  pop      (Intercept) 0.00000 
##  Residual             0.21042 
## Number of obs: 225, groups:  pop:loc, 22; pop, 11
## Fixed Effects:
##    (Intercept)       speciesCL       speciesGA       speciesGR  
##        0.38985        -0.05067         0.03105         0.01338  
##      speciesMR       speciesSC            locO  speciesCL:locO  
##       -0.02840         0.28780         0.15278         0.13277  
## speciesGA:locO  speciesGR:locO  speciesMR:locO  speciesSC:locO  
##       -0.10069         0.23077         0.28893        -0.04707
ggplot(fdat_ind,aes(species,Mycorrhiza,colour=loc))+geom_boxplot()+theme_simple()

ggplot(fdat_ind,aes(pop,Mycorrhiza,colour=loc))+geom_boxplot()+theme_simple()

Mycorrhia model results

  1. Mycorrhiza colonization is significantly higher outside patches of A. petiolata than inside 2a. Species vary in Mycorrhiza colonization, with highest colonization in SC and others relatively similar. 2b. Effect of A. petiolata on Mycorrhiza colonization differs among species, with more pronounced effects for GR and MR, and will little effect for GA
  2. Effect of A. petiolata on Mycorrhiza colonization varies by population, but not as strongly as among species
## Mycorrhiza model (NOTE: No loc*species because not all species present in all populations)
base.Les<-lmer(Lesion~Mycorrhiza*species*loc+(1|pop/loc),data=fdat_ind)
norint.Les<-lmer(Lesion~Mycorrhiza*species*loc+(1|pop),data=fdat_ind)
anova(base.Les,norint.Les)
## refitting model(s) with ML (instead of REML)
norint.Les # inspect coefficients
## Linear mixed model fit by REML ['lmerMod']
## Formula: Lesion ~ Mycorrhiza * species * loc + (1 | pop)
##    Data: fdat_ind
## REML criterion at convergence: -70.9219
## Random effects:
##  Groups   Name        Std.Dev.
##  pop      (Intercept) 0.04934 
##  Residual             0.18543 
## Number of obs: 225, groups:  pop, 11
## Fixed Effects:
##               (Intercept)                 Mycorrhiza  
##                  0.686711                   0.053896  
##                 speciesCL                  speciesGA  
##                 -0.188794                   0.146342  
##                 speciesGR                  speciesMR  
##                  0.090499                  -0.310235  
##                 speciesSC                       locO  
##                 -0.037506                   0.016259  
##      Mycorrhiza:speciesCL       Mycorrhiza:speciesGA  
##                  0.216958                  -0.702730  
##      Mycorrhiza:speciesGR       Mycorrhiza:speciesMR  
##                 -0.020628                   0.559261  
##      Mycorrhiza:speciesSC            Mycorrhiza:locO  
##                  0.136864                  -0.056951  
##            speciesCL:locO             speciesGA:locO  
##                 -0.225004                  -0.389157  
##            speciesGR:locO             speciesMR:locO  
##                 -0.565486                  -0.129568  
##            speciesSC:locO  Mycorrhiza:speciesCL:locO  
##                 -0.007965                   0.098911  
## Mycorrhiza:speciesGA:locO  Mycorrhiza:speciesGR:locO  
##                  0.784046                   0.433494  
## Mycorrhiza:speciesMR:locO  Mycorrhiza:speciesSC:locO  
##                 -0.172551                  -0.051787
Les2<-lmer(Lesion~Mycorrhiza*species*loc-Mycorrhiza:species:loc+(1|pop),data=fdat_ind)
anova(norint.Les,Les2)
## refitting model(s) with ML (instead of REML)
Les3a<-lmer(Lesion~Mycorrhiza*species+Mycorrhiza*loc+(1|pop),data=fdat_ind)
Les3b<-lmer(Lesion~Mycorrhiza*species+species*loc+(1|pop),data=fdat_ind)
Les3c<-lmer(Lesion~species*loc+Mycorrhiza*loc+(1|pop),data=fdat_ind)
anova(norint.Les,Les3a) # TEST species*loc 
## refitting model(s) with ML (instead of REML)
anova(norint.Les,Les3b) # TEST Myc*loc 
## refitting model(s) with ML (instead of REML)
anova(norint.Les,Les3c) # TEST Myc*species
## refitting model(s) with ML (instead of REML)
ggplot(fdat_ind,aes(Mycorrhiza,Lesion))+facet_wrap(~species)+
  geom_point(aes(col=loc))+theme_simple()+geom_smooth(method="lm")

** Lesions model results**

  1. Root lesions increase with mycorrhiza colonization, contrary to predictions
  2. Effect of Mycorrhiza on Lesions depends strongly on species (Myc*species interactions)

Hypothesis 1: Lesion ~ Mycorriza

GIVEN THE NEW ANALYSIS ABOVE, THIS ISN’T WORTHWILE BECAUSE THERE IS NO SIGNIFICANT Myc:species:loc

# calculating the r
indivCODE <- unique(fdata$indiv2)
dfcorr <-data.frame(indivCODE=c(0),r=c(0))
for (i in 1:length(indivCODE)){
  dat1<-NULL
  dat1<-filter(fdata, fdata$indiv2==indivCODE[i])
  dfcorr[i,1] <- indivCODE[i]
  dfcorr[i,2] <- cor(as.numeric(dat1$Mycorrhiza),as.numeric(dat1$Lesion))
}
## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero

## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero

## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero

## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero

## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero

## Warning in cor(as.numeric(dat1$Mycorrhiza), as.numeric(dat1$Lesion)): the
## standard deviation is zero
dfcorr %>% drop_na()->dfcorr

Co-occurence data

#Normal?
car::qqp(dfcorr$r, "norm")
shapiro.test(dfcorr$r)
hist(dfcorr$r)
# Cannot be normal because it is only between -1 and 1 and our mean is at 0.
descdist(dfcorr$r, discrete = FALSE, boot = 500)

# weibull when transformed to be positive?
dfcorr$r.t <- dfcorr$r+1
fit.weibull <- fitdist(dfcorr$r.t, "weibull")
car::qqp(dfcorr$r.t, "weibull", shape = fit.weibull$estimate[[1]])

dfcorr$pop <- gsub("^([0-9]+)[I,O][A-Z]+.*","\\1",dfcorr$indivCODE)
dfcorr$location <- gsub("^[0-9]+([I,O])[A-Z]+.*","\\1",dfcorr$indivCODE)
dfcorr$species <- sub("^[0-9]+[I,O]([A-Z]+).*","\\1",dfcorr$indivCODE)
dfcorr$location <- as.factor(dfcorr$location)
dfcorr$species <- as.factor(dfcorr$species)
dfcorr$pop <- as.factor(dfcorr$pop)

meanR<-mean(dfcorr$r)
SEMR<-sd(dfcorr$r)/sqrt(length(dfcorr$r))
UCI<- meanR+SEMR*1.96
LCI<- meanR-SEMR*1.96
mean(filter(dfcorr, dfcorr$location=="I")$r)
mean(filter(dfcorr, dfcorr$location=="O")$r)


baysian_corrGLMM <- brm(r.t ~ location * species + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
baysian_corrGLMM2 <- brm(r.t ~ location + species + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
baysian_corrGLMM3 <- brm(r.t ~ location + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))
baysian_corrGLMM_NULL <- brm(r.t ~ 1 + (1|pop), data=dfcorr, family = weibull(link="log"), control = list(adapt_delta = 0.9))

waic.BCG_1<-waic(baysian_corrGLMM)
waic.BCG_2<-waic(baysian_corrGLMM2)
waic.BCG_3<-waic(baysian_corrGLMM3)
waic.BCG_NULL<-waic(baysian_corrGLMM_NULL)
compare(waic.BCG_1,waic.BCG_2) #positive means second is better
compare(waic.BCG_2,waic.BCG_3) #negative means first is better
compare(waic.BCG_2,waic.BCG_NULL) #negative means first is better
summary(baysian_corrGLMM2)
bayes_R2(baysian_corrGLMM2)


ggplot(dfcorr, aes(y=r,x=location))+geom_boxplot()+theme_simple()
ggplot(dfcorr, aes(y=r,x=species))+geom_boxplot()+theme_simple()

Totals data

# model 
GLMMmyc_lesion <- glmer(binlesion ~ None.Myc * species + (1 | pop) , data = ndatas,family = binomial(link = "logit"))
summary(GLMMmyc_lesion)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: binlesion ~ None.Myc * species + (1 | pop)
##    Data: ndatas
## 
##      AIC      BIC   logLik deviance df.resid 
##   5388.9   5433.4  -2681.5   5362.9      213 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.9159 -3.2800 -0.7426  2.8805 11.6817 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  pop    (Intercept) 0.07038  0.2653  
## Number of obs: 226, groups:  pop, 11
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -0.750276   0.090719  -8.270  < 2e-16 ***
## None.Myc           -0.017781   0.038118  -0.466 0.640866    
## speciesCL           0.415585   0.056183   7.397 1.39e-13 ***
## speciesGA           0.444280   0.053486   8.306  < 2e-16 ***
## speciesGR           0.266576   0.073678   3.618 0.000297 ***
## speciesMR           0.146355   0.057804   2.532 0.011344 *  
## speciesSC           0.007849   0.051574   0.152 0.879040    
## None.Myc:speciesCL  0.156587   0.048414   3.234 0.001219 ** 
## None.Myc:speciesGA  0.191932   0.050731   3.783 0.000155 ***
## None.Myc:speciesGR  0.338206   0.069696   4.853 1.22e-06 ***
## None.Myc:speciesMR  0.063131   0.053541   1.179 0.238360    
## None.Myc:speciesSC  0.210590   0.051242   4.110 3.96e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Nn.Myc spcsCL spcsGA spcsGR spcsMR spcsSC N.M:CL N.M:GA
## None.Myc    -0.134                                                        
## speciesCL   -0.369  0.213                                                 
## speciesGA   -0.385  0.223  0.649                                          
## speciesGR   -0.279  0.160  0.466  0.500                                   
## speciesMR   -0.340  0.207  0.572  0.614  0.438                            
## speciesSC   -0.381  0.233  0.627  0.660  0.477  0.588                     
## Nn.Myc:spCL  0.098 -0.786 -0.235 -0.155 -0.110 -0.155 -0.171              
## Nn.Myc:spGA  0.094 -0.750 -0.141 -0.235 -0.116 -0.146 -0.169  0.593       
## Nn.Myc:spGR  0.086 -0.551 -0.142 -0.151 -0.111 -0.128 -0.145  0.431  0.406
## Nn.Myc:spMR  0.090 -0.711 -0.144 -0.148 -0.092 -0.111 -0.159  0.560  0.534
## Nn.Myc:spSC  0.116 -0.750 -0.191 -0.195 -0.145 -0.169 -0.045  0.589  0.555
##             N.M:GR N.M:MR
## None.Myc                 
## speciesCL                
## speciesGA                
## speciesGR                
## speciesMR                
## speciesSC                
## Nn.Myc:spCL              
## Nn.Myc:spGA              
## Nn.Myc:spGR              
## Nn.Myc:spMR  0.394       
## Nn.Myc:spSC  0.429  0.527
options(na.action="na.fail")
dredge(GLMMmyc_lesion)
## Fixed term is "(Intercept)"
r.squaredGLMM(GLMMmyc_lesion)
## The result is correct only if all data used by the model has not changed since model was fitted.
##        R2m        R2c 
## 0.02739802 0.04008648
# will need to change this ASAP

ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path))+
  stat_summary(fun.y=mean, geom="point")+
  geom_smooth(method = "lm")+
  geom_point()+
  theme_simple()+
  facet_wrap("species")+
  labs(x = "Mycorrhizal colonization", y="Lesions")

Hypothesis 2: Mycorrhiza & Lesion ~ invasion

Does location impact the effect of myc on lesion?

GLMMmyc_lesion_loc <- glmer(binlesion ~ None.Myc * location * species + (1 | pop), data = ndatas, family = binomial(link = "logit"))
summary(GLMMmyc_lesion_loc)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: binlesion ~ None.Myc * location * species + (1 | pop)
##    Data: ndatas
## 
##      AIC      BIC   logLik deviance df.resid 
##   5185.5   5271.0  -2567.7   5135.5      201 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.6366 -3.2043 -0.3462  2.6760 13.6450 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  pop    (Intercept) 0.0698   0.2642  
## Number of obs: 226, groups:  pop, 11
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -0.725383   0.100890  -7.190 6.49e-13 ***
## None.Myc                     -0.062764   0.054987  -1.141 0.253688    
## locationO                    -0.068081   0.081316  -0.837 0.402458    
## speciesCL                     0.200476   0.081638   2.456 0.014063 *  
## speciesGA                     0.296660   0.077883   3.809 0.000139 ***
## speciesGR                    -1.181296   0.193649  -6.100 1.06e-09 ***
## speciesMR                     0.277361   0.084158   3.296 0.000982 ***
## speciesSC                    -0.021562   0.073106  -0.295 0.768042    
## None.Myc:locationO            0.093031   0.080266   1.159 0.246440    
## None.Myc:speciesCL            0.324299   0.071788   4.517 6.26e-06 ***
## None.Myc:speciesGA            0.189403   0.072737   2.604 0.009216 ** 
## None.Myc:speciesGR           -0.300771   0.187873  -1.601 0.109393    
## None.Myc:speciesMR            0.255496   0.077895   3.280 0.001038 ** 
## None.Myc:speciesSC            0.101743   0.070113   1.451 0.146743    
## locationO:speciesCL           0.398803   0.106396   3.748 0.000178 ***
## locationO:speciesGA           0.301929   0.101835   2.965 0.003028 ** 
## locationO:speciesGR           1.889287   0.212942   8.872  < 2e-16 ***
## locationO:speciesMR          -0.260917   0.111010  -2.350 0.018754 *  
## locationO:speciesSC           0.143608   0.103149   1.392 0.163850    
## None.Myc:locationO:speciesCL -0.282839   0.100392  -2.817 0.004842 ** 
## None.Myc:locationO:speciesGA -0.003347   0.105512  -0.032 0.974693    
## None.Myc:locationO:speciesGR  0.463051   0.207868   2.228 0.025906 *  
## None.Myc:locationO:speciesMR -0.384191   0.109499  -3.509 0.000450 ***
## None.Myc:locationO:speciesSC  0.233745   0.102953   2.270 0.023183 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
emm<-emmeans(GLMMmyc_lesion_loc, list(~None.Myc| species | location, ~None.Myc| location), type = "response" )
## NOTE: Results may be misleading due to involvement in interactions
summary(emm)
## $`emmeans of None.Myc, location | species`
## species = AH:
##      None.Myc location      prob         SE df  asymp.LCL asymp.UCL
##  1.058488e-16 I        0.3262087 0.02217532 NA 0.28432176 0.3710671
##  1.058488e-16 O        0.3114255 0.02097601 NA 0.27186461 0.3539447
## 
## species = CL:
##      None.Myc location      prob         SE df  asymp.LCL asymp.UCL
##  1.058488e-16 I        0.3717055 0.02222786 NA 0.32927810 0.4162074
##  1.058488e-16 O        0.4516058 0.02288001 NA 0.40727139 0.4967219
## 
## species = GA:
##      None.Myc location      prob         SE df  asymp.LCL asymp.UCL
##  1.058488e-16 I        0.3944313 0.02175155 NA 0.35269620 0.4377651
##  1.058488e-16 O        0.4514349 0.02239380 NA 0.40803107 0.4955905
## 
## species = GR:
##      None.Myc location      prob         SE df  asymp.LCL asymp.UCL
##  1.058488e-16 I        0.1293544 0.02255808 NA 0.09118441 0.1803321
##  1.058488e-16 O        0.4786449 0.02657122 NA 0.42698389 0.5307668
## 
## species = MR:
##      None.Myc location      prob         SE df  asymp.LCL asymp.UCL
##  1.058488e-16 I        0.3898311 0.02331380 NA 0.34521822 0.4363672
##  1.058488e-16 O        0.3149625 0.02075782 NA 0.27576180 0.3569890
## 
## species = SC:
##      None.Myc location      prob         SE df  asymp.LCL asymp.UCL
##  1.058488e-16 I        0.3214874 0.01942629 NA 0.28465376 0.3606841
##  1.058488e-16 O        0.3381795 0.02091737 NA 0.29847025 0.3803078
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $`emmeans of None.Myc | location`
## location = I:
##      None.Myc      prob         SE df asymp.LCL asymp.UCL
##  1.058488e-16 0.3107155 0.01873775 NA 0.2752289 0.3485770
## 
## location = O:
##      None.Myc      prob         SE df asymp.LCL asymp.UCL
##  1.058488e-16 0.3887098 0.01966447 NA 0.3509329 0.4278726
## 
## Results are averaged over the levels of: species 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale
ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path))+
  geom_smooth(method = "lm")+
  geom_point(aes(color=location))+
  theme_simple()+
  facet_wrap("species")+
  labs(x = "Mycorrhizal colonization", y="Lesions")

pca <- prcomp(ndata[,c(5,9)])
location <- ndata[,13]
g <- ggbiplot(pca, obs.scale = 1, var.scale = 1, 
                groups = location, ellipse = TRUE, 
                circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal', 
                 legend.position = 'top')
print(g)

Does location impact lesion alone?

glmerLESION_disp <- glmer (binlesion ~ location * species + (1|pop) + (1|indiv), data = ndatas, family = binomial(link="logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00135513 (tol =
## 0.001, component 1)
summary(glmerLESION_disp)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: binlesion ~ location * species + (1 | pop) + (1 | indiv)
##    Data: ndatas
## 
##      AIC      BIC   logLik deviance df.resid 
##   2011.7   2059.6   -991.8   1983.7      212 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.53231 -0.13210 -0.00971  0.11216  1.14859 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  indiv  (Intercept) 1.11959  1.0581  
##  pop    (Intercept) 0.06426  0.2535  
## Number of obs: 226, groups:  indiv, 226; pop, 11
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)   
## (Intercept)         -0.85482    0.30710  -2.784  0.00538 **
## locationO           -0.10392    0.39475  -0.263  0.79235   
## speciesCL            0.26989    0.39300   0.687  0.49225   
## speciesGA            0.38852    0.37584   1.034  0.30127   
## speciesGR           -0.88023    0.62967  -1.398  0.16214   
## speciesMR            0.38516    0.41882   0.920  0.35776   
## speciesSC           -0.09661    0.35549  -0.272  0.78581   
## locationO:speciesCL  0.45477    0.52759   0.862  0.38870   
## locationO:speciesGA  0.47505    0.50259   0.945  0.34456   
## locationO:speciesGR  1.77188    0.77503   2.286  0.02224 * 
## locationO:speciesMR -0.28686    0.55694  -0.515  0.60651   
## locationO:speciesSC -0.02068    0.48262  -0.043  0.96582   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) loctnO spcsCL spcsGA spcsGR spcsMR spcsSC lcO:CL lcO:GA
## locationO   -0.703                                                        
## speciesCL   -0.741  0.551                                                 
## speciesGA   -0.778  0.579  0.616                                          
## speciesGR   -0.455  0.347  0.354  0.376                                   
## speciesMR   -0.682  0.522  0.539  0.570  0.340                            
## speciesSC   -0.807  0.608  0.637  0.669  0.393  0.588                     
## lctnO:spcCL  0.524 -0.747 -0.716 -0.431 -0.253 -0.390 -0.453              
## lctnO:spcGA  0.566 -0.787 -0.446 -0.727 -0.272 -0.411 -0.488  0.587       
## lctnO:spcGR  0.354 -0.513 -0.273 -0.288 -0.807 -0.269 -0.307  0.378  0.400
## lctnO:spcMR  0.502 -0.709 -0.396 -0.411 -0.250 -0.732 -0.434  0.530  0.558
## lctnO:spcSC  0.575 -0.818 -0.450 -0.473 -0.283 -0.428 -0.721  0.612  0.643
##             lcO:GR lcO:MR
## locationO                
## speciesCL                
## speciesGA                
## speciesGR                
## speciesMR                
## speciesSC                
## lctnO:spcCL              
## lctnO:spcGA              
## lctnO:spcGR              
## lctnO:spcMR  0.366       
## lctnO:spcSC  0.419  0.580
## convergence code: 0
## Model failed to converge with max|grad| = 0.00135513 (tol = 0.001, component 1)
options(na.action="na.fail")
dredge(glmerLESION_disp)
## Fixed term is "(Intercept)"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00135513 (tol =
## 0.001, component 1)
r.squaredGLMM(glmerLESION_disp)
##        R2m        R2c 
## 0.03076813 0.04469034
ggplot(aes(x = location, y = 100-None.Path, color=location), data=ndata)+
  geom_boxplot()+
  geom_jitter(width=0.1)+
  facet_wrap("species")+
  theme_simple()

Does location impact mycorrhiza alone?

glmerMYC_disp <- glmer (binmycorr ~ location * species + (1|pop) + (1|indiv), data = ndata, family = binomial(link="logit"))
summary(glmerMYC_disp)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: binmycorr ~ location * species + (1 | pop) + (1 | indiv)
##    Data: ndata
## 
##      AIC      BIC   logLik deviance df.resid 
##   2073.3   2121.2  -1022.7   2045.3      212 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.35863 -0.12130  0.01523  0.09795  1.40036 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  indiv  (Intercept) 2.062e+00 1.4361234
##  pop    (Intercept) 1.044e-08 0.0001022
## Number of obs: 226, groups:  indiv, 226; pop, 11
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)          0.121220   0.390387   0.310   0.7562  
## locationO           -0.008745   0.526879  -0.017   0.9868  
## speciesCL           -0.200766   0.515746  -0.389   0.6971  
## speciesGA           -0.375379   0.491194  -0.764   0.4447  
## speciesGR           -1.324367   0.829848  -1.596   0.1105  
## speciesMR           -0.505583   0.551771  -0.916   0.3595  
## speciesSC           -0.954255   0.471189  -2.025   0.0428 *
## locationO:speciesCL  0.053276   0.706272   0.075   0.9399  
## locationO:speciesGA  0.233201   0.669292   0.348   0.7275  
## locationO:speciesGR  0.990565   1.025133   0.966   0.3339  
## locationO:speciesMR -0.245112   0.745205  -0.329   0.7422  
## locationO:speciesSC -0.457001   0.645298  -0.708   0.4788  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) loctnO spcsCL spcsGA spcsGR spcsMR spcsSC lcO:CL lcO:GA
## locationO   -0.741                                                        
## speciesCL   -0.757  0.561                                                 
## speciesGA   -0.795  0.589  0.602                                          
## speciesGR   -0.470  0.349  0.356  0.374                                   
## speciesMR   -0.708  0.524  0.536  0.562  0.333                            
## speciesSC   -0.829  0.614  0.627  0.659  0.390  0.586                     
## lctnO:spcCL  0.553 -0.746 -0.730 -0.439 -0.260 -0.391 -0.458              
## lctnO:spcGA  0.583 -0.787 -0.442 -0.734 -0.274 -0.413 -0.483  0.587       
## lctnO:spcGR  0.381 -0.514 -0.288 -0.303 -0.809 -0.269 -0.316  0.383  0.405
## lctnO:spcMR  0.524 -0.707 -0.396 -0.416 -0.246 -0.740 -0.434  0.527  0.557
## lctnO:spcSC  0.605 -0.816 -0.458 -0.481 -0.285 -0.428 -0.730  0.609  0.643
##             lcO:GR lcO:MR
## locationO                
## speciesCL                
## speciesGA                
## speciesGR                
## speciesMR                
## speciesSC                
## lctnO:spcCL              
## lctnO:spcGA              
## lctnO:spcGR              
## lctnO:spcMR  0.363       
## lctnO:spcSC  0.420  0.577
options(na.action="na.fail")
dredge(glmerMYC_disp)
## Fixed term is "(Intercept)"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0029301 (tol =
## 0.001, component 1)
r.squaredGLMM(glmerMYC_disp)
##       R2m       R2c 
## 0.0410656 0.0410656
ggplot(aes(x = location, y = 100-None.Myc, color=location), data=ndata)+
  geom_boxplot()+
  geom_jitter(width=0.1)+
  facet_wrap("species")+
  theme_simple()

Hypothesis 3 + 4: Lesion ~ Mycorriza * invasion

ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path, color=location))+
  geom_smooth(method = "lm")+
  geom_point()+
  theme_simple()+
  facet_wrap("species")+
  labs(x = "Mycorrhizal colonization", y="Lesions")

ggplot(data = ndata,aes(x = 1-None.Myc, y = 100-None.Path, color=location))+
  geom_smooth(method = "lm")+
  geom_point()+
  theme_simple()+
  labs(x = "Mycorrhizal colonization", y="Lesions")

Bring things together

Soil_characteristics$Location<-as.character(Soil_characteristics$Location)
for(i in 1:length(Soil_characteristics$Location))if(Soil_characteristics[i,2]=="in"){
   Soil_characteristics[i,2]<-"I"
}else{
  Soil_characteristics[i,2]<-"O"
}
Soil_characteristics$Location<-as.factor(Soil_characteristics$Location)  
names(Soil_characteristics)[1:2]<-c("pop","location")

fulldata <- merge(Soil_characteristics, ndata, by=c("pop","location"))
# note i used the SCALED dataframe for the total histology
fulldata <- merge(fulldata, Root_size, by=c("indiv","location","pop","species"))
for(i in 1:length(fulldata$indiv)){
fulldata$mean_rootsize[i]<-(sum(fulldata[i,20:29])/length(fulldata[,20:29]))
}
fulldata<-fulldata[,-c(20:29)]

binlesionF<-cbind(fulldata$None.Path,100-fulldata$None.Path)
binmycorrF<-cbind(fulldata$None.Myc,100-fulldata$None.Myc)
binherbivoryF<-cbind(fulldata$None.Herb,100-fulldata$None.Herb)


fulldata$pop<-as.factor(fulldata$pop)
# model myc_lesion 
mod1 <- glmmadmb(binlesionF ~ scale(None.Myc)  * scale(mean_rootsize) * species + scale(pH) + scale(Agg_stability) + scale(N_percent) + scale(C_percent) + (1 | pop), data = fulldata,family = "binomial")
summary(mod1)
## 
## Call:
## glmmadmb(formula = binlesionF ~ scale(None.Myc) * scale(mean_rootsize) * 
##     species + scale(pH) + scale(Agg_stability) + scale(N_percent) + 
##     scale(C_percent) + (1 | pop), data = fulldata, family = "binomial")
## 
## AIC: 3914.8 
## 
## Coefficients:
##                                                Estimate Std. Error z value
## (Intercept)                                     -0.9247     0.1360   -6.80
## scale(None.Myc)                                 -0.0167     0.0601   -0.28
## scale(mean_rootsize)                             0.1370     0.0727    1.88
## speciesCL                                        0.4485     0.0730    6.15
## speciesGA                                        1.7767     0.1557   11.41
## speciesGR                                        0.3890     0.1028    3.78
## speciesMR                                        0.3305     0.0742    4.45
## speciesSC                                        0.4228     0.0721    5.87
## scale(pH)                                       -0.1158     0.0262   -4.43
## scale(Agg_stability)                            -0.0166     0.0322   -0.51
## scale(N_percent)                                -0.6563     0.1253   -5.24
## scale(C_percent)                                 0.6935     0.1489    4.66
## scale(None.Myc):scale(mean_rootsize)             0.1930     0.0619    3.12
## scale(None.Myc):speciesCL                        0.0459     0.0712    0.64
## scale(None.Myc):speciesGA                        0.3262     0.1761    1.85
## scale(None.Myc):speciesGR                        0.0850     0.1080    0.79
## scale(None.Myc):speciesMR                        0.1156     0.0755    1.53
## scale(None.Myc):speciesSC                        0.3634     0.0793    4.58
## scale(mean_rootsize):speciesCL                   0.0602     0.0806    0.75
## scale(mean_rootsize):speciesGA                   1.0349     0.1636    6.33
## scale(mean_rootsize):speciesGR                  -0.3641     0.1051   -3.46
## scale(mean_rootsize):speciesMR                  -0.0544     0.0882   -0.62
## scale(mean_rootsize):speciesSC                   0.0131     0.1275    0.10
## scale(None.Myc):scale(mean_rootsize):speciesCL  -0.0199     0.0656   -0.30
## scale(None.Myc):scale(mean_rootsize):speciesGA  -0.1636     0.1831   -0.89
## scale(None.Myc):scale(mean_rootsize):speciesGR  -0.4582     0.1258   -3.64
## scale(None.Myc):scale(mean_rootsize):speciesMR  -0.2581     0.0762   -3.39
## scale(None.Myc):scale(mean_rootsize):speciesSC  -0.1735     0.1278   -1.36
##                                                Pr(>|z|)    
## (Intercept)                                     1.1e-11 ***
## scale(None.Myc)                                 0.78093    
## scale(mean_rootsize)                            0.05950 .  
## speciesCL                                       7.9e-10 ***
## speciesGA                                       < 2e-16 ***
## speciesGR                                       0.00016 ***
## speciesMR                                       8.4e-06 ***
## speciesSC                                       4.5e-09 ***
## scale(pH)                                       9.6e-06 ***
## scale(Agg_stability)                            0.60758    
## scale(N_percent)                                1.6e-07 ***
## scale(C_percent)                                3.2e-06 ***
## scale(None.Myc):scale(mean_rootsize)            0.00181 ** 
## scale(None.Myc):speciesCL                       0.51893    
## scale(None.Myc):speciesGA                       0.06396 .  
## scale(None.Myc):speciesGR                       0.43096    
## scale(None.Myc):speciesMR                       0.12571    
## scale(None.Myc):speciesSC                       4.5e-06 ***
## scale(mean_rootsize):speciesCL                  0.45510    
## scale(mean_rootsize):speciesGA                  2.5e-10 ***
## scale(mean_rootsize):speciesGR                  0.00053 ***
## scale(mean_rootsize):speciesMR                  0.53735    
## scale(mean_rootsize):speciesSC                  0.91845    
## scale(None.Myc):scale(mean_rootsize):speciesCL  0.76130    
## scale(None.Myc):scale(mean_rootsize):speciesGA  0.37155    
## scale(None.Myc):scale(mean_rootsize):speciesGR  0.00027 ***
## scale(None.Myc):scale(mean_rootsize):speciesMR  0.00071 ***
## scale(None.Myc):scale(mean_rootsize):speciesSC  0.17451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of observations: total=189, pop=9 
## Random effect variance(s):
## Warning in .local(x, sigma, ...): 'sigma' and 'rdig' arguments are present
## for compatibility only: ignored
## Group=pop
##             Variance StdDev
## (Intercept)   0.1348 0.3671
## 
## 
## Log-likelihood: -1928.4
mod1.dredge <- dredge(mod1)
## Fixed term is "(Intercept)"
mod1.dredgegood <- subset(mod1.dredge,mod1.dredge$delta<=7)
# model myc 
mod2 <- glmmadmb(binmycorrF ~ scale(mean_rootsize) * species + scale(pH) + scale(Agg_stability) + scale(N_percent) + scale(C_percent) + (1 | pop) , data = fulldata,family = "binomial")
summary(mod2)
## 
## Call:
## glmmadmb(formula = binmycorrF ~ scale(mean_rootsize) * species + 
##     scale(pH) + scale(Agg_stability) + scale(N_percent) + scale(C_percent) + 
##     (1 | pop), data = fulldata, family = "binomial")
## 
## AIC: 6440.8 
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    -0.00125    0.10116   -0.01  0.99014    
## scale(mean_rootsize)            0.12660    0.06462    1.96  0.05008 .  
## speciesCL                      -0.04954    0.06748   -0.73  0.46286    
## speciesGA                      -0.03748    0.13458   -0.28  0.78065    
## speciesGR                      -0.40656    0.08369   -4.86  1.2e-06 ***
## speciesMR                      -0.25895    0.06867   -3.77  0.00016 ***
## speciesSC                      -0.82100    0.06393  -12.84  < 2e-16 ***
## scale(pH)                       0.08393    0.02423    3.46  0.00053 ***
## scale(Agg_stability)            0.18425    0.03032    6.08  1.2e-09 ***
## scale(N_percent)                0.14058    0.10894    1.29  0.19691    
## scale(C_percent)               -0.14612    0.13142   -1.11  0.26619    
## scale(mean_rootsize):speciesCL -0.17833    0.07173   -2.49  0.01292 *  
## scale(mean_rootsize):speciesGA -0.17024    0.14277   -1.19  0.23309    
## scale(mean_rootsize):speciesGR -0.26173    0.08015   -3.27  0.00109 ** 
## scale(mean_rootsize):speciesMR -0.12646    0.07473   -1.69  0.09062 .  
## scale(mean_rootsize):speciesSC -1.26599    0.12042  -10.51  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of observations: total=189, pop=9 
## Random effect variance(s):
## Warning in .local(x, sigma, ...): 'sigma' and 'rdig' arguments are present
## for compatibility only: ignored
## Group=pop
##             Variance StdDev
## (Intercept)   0.0652 0.2554
## 
## 
## Log-likelihood: -3203.42
mod2.dredge <- dredge(mod2)
## Fixed term is "(Intercept)"
mod2.dredgegood <- subset(mod2.dredge,mod2.dredge$delta<=7)

Molecular data Analysis

Hypothesis 1: Mycorrhiza_div ~ invasion

Hypothesis 2: Mycorrhiza_div ~ Pathoge_div